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Chapter  1  Introduction 

1.1  Problem  Description 

The  purpose  of  this  report  is  to  develop  automatic  aircraft  flight  control  systems  using 
the  concept  of  Total  Energy  Control  System  (TECS)  in  the  lateral  directional  axis  for 
the  McDonnell  Douglas  F- 15  Eagle.  Nonlinear  simulation  was  originally  developed 
from  actual  flight  tests  conducted  by  the  U.S.  Air  Force  and  generated  into  a  Fortran 
code.  Lt.  James  P.  Dutton  then  took  the  flight  data  code  and  generated  a  SIMULINK 
file  [10].  Using  the  SIMULINK  file  a  control  law  developed  on  the  linearized  model 
can  then  be  tested  on  the  nonlinear  model. 

1.2  Analysis  of  Data  Provided 

Due  to  the  complexity  of  the  flight  data  developed  for  the  F-15,  the  next  sections  are 
used  to  describe  important  support  modules  contained  in  [5].  The  numerical  data 
were  provided  in  the  format  of  a  Genesis  simulation  used  at  Wright-Patterson  AFB 
and  the  coding  was  subsequently  converted  from  FORTRAN  to  MATLAB.  A  de¬ 
tailed  description  of  the  model  is  located  in  Chapter  2. 

Evaluation  of  the  open-loop  nonlinear  model  and  design  of  the  TECS  control  law 
will  be  accomplished  at  two  distinct  flight  points  shown  in  Table  1.1.  Other  flight 
conditions  throughout  the  flight  envelope  can  be  investigated  in  future  study. 

Table  1.1:  Selected  Flight  Conditions 


Flight 

Altitude 

Vtas 

Mach 

Point 

(ft) 

(ft/s) 

Number 

1 

9,800 

539.1 

0.5 

2 

30,000 

497.3 

0.5 

1.2.1  Model  Characteristics 

The  model  is  an  integration  of  several  modules,  each  performing  a  specific  function. 
Included  in  these  modules  are  the  aerodynamic  data,  the  propulsion  data,  the  atmo¬ 
spheric  conditions,  as  well  as  the  nonlinear  equations  of  motion  for  the  F-15.  As  far 
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as  the  modeling  of  the  control  actuators  and  sensors  are  concerned,  they  can  be  en¬ 
tered  into  the  model  using  the  capabilities  within  SIMULINK  to  create  linear  filter 
dynamics.  Figure  1.1  summarizes  the  integration  of  the  modules  that  form  the  com¬ 
plete  system.  -  . 


Figure  1.1  Integration  of  System  Module  Components 

The  aircraft  is  the  McDonnell  F-15  Eagle,  the  first  line  of  defense  for  operational 
high-performance  aircraft.  It  is  powered  by  two  afterburning  turbofan  engines,  each 
providing  approximately  32,000  pounds  of  thrust  at  maximum  power.  The  primary 
flight  controls  are  the  horizontal  stabilizers  capable  of  both  symmetric  and  differen¬ 
tial  movement,  conventional  ailerons,  and  dual  vertical  rudders.  There  are  a  total  of 
six  actuators,  one  for  each  surface  mentioned  above.  All  actuators  are  modeled  iden¬ 
tically  with  rate  limits  of  24  deg/sec  and  first-order  response  characteristics  of 

^ 

Position  limitations  on  the  actual  component  and  the  sign  conventions  for  posi¬ 
tive  deflections  are  summarized  in  Table  1.2.  The  aircraft  mass  and  geometry 


characteristics  are  summarized  in  Table  1.3. 


Table  1.2:  Flight  Control  Surfaces  Characteristics 


Control 

Surface 

Symbol 

Limits 

Sign 

Convention 

(+) 

Symmetric 

Stabilator 

Sh 

t20° 

Trailing  edge 
down 

Differential 

Stabilator 

5d 

l5°-25° 

Left  trailing 
edge  down 

Aileron 

5a 

±20° 

Left  trailing 
edge  down 

Sr 


Rudder 


±30° 


Trailing  edge 
left 


Table  1.3:  F-15  Mass  and  Geometry  Characteristics 


Parameter 

Symbol 

IB^I 

Value 

Wing  Area 

S 

_^Wing  span 

b 

ft 

42.8 

Mean  aero¬ 
dynamic 
chord 

c 

ft 

15.95 

Aircraft 

Weight 

W 

lb 

45,000 

Moments 

Ix 

slug-ft^ 

28,700 

of 

h 

slug-ft^ 

165,100 

Inertia 

Iz 

slug-ft^ 

187,900 

Products 

^xz 

slug-ft^ 

-520 

of 

^xy 

slug-ft^ 

0 

Inertia 

^yz 

slug-ft^ 

0 

1.2.2  Aerodynamic  Model 

The  aerodynamic  model  is  composed  of  multidimensional  table  with  interpolation  to 
form  nonlinear  function  generators.  By  using  a  vast  database  the  highly  nonlinear 
aerodynamics  encountered  in  the  large  flight  envelope  of  the  F-15  are  properly  rep¬ 
resented.  This  will  be  important  for  evaluating  the  TECS  controller  being  developed 
for  the  linear  model  on  the  nonlinear  system.  Most  of  the  aerodynamic  qualities  are 
determined  from  the  Mach  number.  In  addition,  combination  of  the  sideslip  angle  p, 
angle  of  attack  a,  and  the  symmetric  stabilator  deflection  5^  are  also  used  to  deter¬ 
mine  the  flight  environment. 

The  task  of  the  aerodynamic  module  is  to  provide  the  nondimensionalized  force 
and  moment  coefficients,  which  are  then  used  to  calculate  the  associated  moments 
and  forces  of  the  F-15.  The  basic  equations  used  for  the  coefficients  are 


Coefficients  of  forces: 


~  ^Lo'^  ^^Lnz^z 
~  ^Do  ^^Dalt  ^^Dnoz 
CyiTr  ^Yo  ^Y&A^A  ^Y&D^D  ~  ^^YSR^R 
Coefficients  of  moments: 

-  ^lo  +  ^IbA^A  ■*■  ^ISD^D  ~  ^^IBR^SR  ^  ^  ^IpP  ^/r^) 

~  ^mo  ^^mnz^z  ^  ^  ^mq9  + 

The  terms  C,  AC,  AN,  and  K  are  outputs  from  the  function  generation  routines,  and 
are  either  read  directly  or  are  derived  from  linear  interpolation  of  tabular  data.  The 
parameters  affecting  each  of  the  coefficients  are  located  in  Table  1.4. 
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Table  1.4:  Independent  Parameters  Affecting  the  Aerodynamic  Coefficients 


Aero 

Coefficients 

Independent 

Parameters 

Aero 

Coefficients 

Independent 

Parameters 

o 

1-1 

U 

M,  a,  5h 

Clp 

M,  a 

^^Lnz 

M 

Cir 

M,  a 

c 

'^mo 

M,  a,  Sjj 

^no 

M,  a,  p 

^^mnz 

M 

^n8a 

M,  a 

^mq 

M,a 

Cn8D 

M,  a 

^ma 

M,  a 

^Qi8r 

M,  a,  p 

ANo 

M 

K^8Rn 

M,  a 

Cd 

KgRl 

M 

(a<32) 

Clo*  M 

^np 

M,  a 

(32<a<40) 

Cio,  M,  a 

Cm 

M,  a 

(a  >40) 

^Lo’  ® 

^yo 

M.  a,p 

h 

^y8A 

M,  a 

^^Dnoz 

M,  5pLA 

^y8D 

M,  a 

Cio 

M,a,P 

^^y8R 

M,  a,5R 

M,a,5R 

K^8Ry 

M 

The  total  forces  and  moments  are  calculated  from  the  following  equations 
L  =  qSCi^ 

D  =  qSCj) 

Y  =  qSCy 
and 

ZL=  qSbCj 
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IM=qScC^ 

m=qSbC„ 

where  q=l/2  pV^  is  the  dynamic  pressure,  S  is  the  wing  planform  area,  b  is  the  wing 
span,  and  c  is  the  mean  aerodynamic  chord. 

1.2.3  Propulsion  Model 

The  propulsion  module  has  two  distinct  sub-modules,  one  for  each  engine.  Although 
each  engine  is  similar,  they  will  not  produce  the  same  amount  of  thrust  for  the  same 
power  settings,  thus  the  need  for  two  different  sub-modules.  The  thrust  vectors  are 
aligned  with  the  aircraft’s  x-axis.  The  thrust  produced  by  each  engine  is  a  function  of 
altitude  h,  Mach  number  M,  and  throtde  setting  5pla-  Furthermore,  each  module  is 
subdivided  into  two  sections-the  core  and  the  afterburner  section,  each  with  its  indi¬ 
vidual  sequencing  logic. 

The  throttle  inputs  to  the  engine  module  are  in  degrees  of  power-level-angle 
(PL A)  with  a  minimum  angle  of  20®  and  a  maximum  position  of  127®.  The  core  sec¬ 
tion  responds  to  the  setting  up  to  83®,  while  the  afterburner  initiates  at  91®.  The  core 
section  has  first-order  dynamics  and  a  rate  limiter  added  in  order  to  simulate  the  spool 
up  time.  The  afterburner  has  a  rate  limiter  and  a  sequencing  logic  to  model  the  fuel 
pump  and  pressure  regulator  effects. 

1.2.4  Atmospheric  Module 

The  atmospheric  model’s  data  is  based  on  tables  from  the  U.S.  Standard  Atmosphere 
(1962).  This  model  calculates  values  for  the  speed  of  sound,  the  acceleration  due  to 
gravity,  air  density,  viscosity,  and  ambient  static  pressure  and  temperature  based  on 
the  aircraft  altitude.  Linear  interpolation  is  used  between  table  values  for  altitudes 
from  0  to  90  km. 

1.2.5  Equations  of  Motion 

The  nonlinear  equations  of  motion  used  in  the  system  model  are  based  on  the  deriva¬ 
tions  by  Duke,  Antoniewicz,  and  Krambeer  in  [6].  These  equations  model  the  six-de- 
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gree-of-freedom  dynamics  of  a  rigid  aircraft  flying  over  a  flat,  non-rotating  Earth. 
These  derivations  are  detailed  in  the  next  chapter. 


Chapter  2  The  Nonlinear  F-15  Model 

Linearized  aircraft  models  are  useful  in  the  early  stage  of  control  law  develop¬ 
ment.  A  design  flight  envelope  can  be  defined  based  on  a  number  of  flight  points. 
However,  the  control  law  will  eventually  have  to  be  validated  in  nonlinear  simula¬ 
tion,  prior  to  the  actual  flight  tests. 

2.1  Derivation  of  Nonlinear  State  Equations 

In  practice,  most  control  laws  are  formulated  by  decoupling  the  longitudinal  motion 
from  the  lateral  motion.  This  is  not  a  correct  practice  in  all  cases.  Certain  simplifying 
assumptions  must  first  be  defined.  The  first  of  which,  the  aircraft  must  be  symmetri¬ 
cal.  In  a  generic  sense,  most  aircraft  are  symmetrical,  however  a  refueling  boom  on 
one  side  of  the  aircraft  such  as  the  A-4  or  a  different  payload  strapped  under  each 
wing  on  the  aircraft  can  cause  this  assumption  to  be  invalid.  The  other  assumption  is 
to  have  the  aircraft  trimmed  at  a  specific  reference  trajectory;  for  instance  straight  and 
level  flight  With  this  specific  reference  condition  and  a  symmetrical  aircraft,  the  lon¬ 
gitudinal  and  lateral  modes  can  always  be  decoupled.  In  the  next  section,  the  aircraft 
will  not  have  either  of  these  assumptions.  Instead,  the  model  will  be  considered  as  a 
rigid  aircraft  of  constant  mass  flying  over  a  flat,  non-rotating  earth. 

2.1.1  Reference  Systems 

There  are  three  primary  reference  systems  associated  with  an  aircraft  in  flight.  These 
are  the  body-axis,  wind-axis,  the  vehicle-carried  vertical-axis  systems.  One  axis  sys¬ 
tem  is  preferred  over  another  axis  depending  on  the  area  of  applications. 

The  rotational  equations  of  motion  are  most  easily  referenced  to  the  body-axis 
system.  The  body-axis  rotational  rates  are  measurable  within  the  aircraft  by  sensors 
fixed  in  the  body  frame.  The  body  axis-system  has  its  origin  at  the  aircraft’s  center 
of  gravity,  the  x-axis  extends  out  the  nose  of  the  aircraft,  the  y-axis  exits  the  right 
wing  of  the  aircraft,  and  the  z-axis  drops  out  the  bottom  of  the  aircraft.  The  positive 
rotation  rates  follow  the  right-hand  rule.  The  x-,  y-,  and  z-  components  of  the  aircraft 
velocity  are  u,  v,  w  respectively.  In  addition  the  roll  rates  are  p,  q,  and  r  with  the  mo- 


merits  being  L,  M,  and  N  for  the  x-,  y-,  and  z-axis  respectively.  Figure  2. 1  shows  these 
relationships. 


Figure  2.1  Aircraft  Body- Axis  System 

The  second-axis  system  is  the  wind-axis  system  which  is  mainly  used  for  the 
translational  equations  of  motion.  Due  to  the  fact  that  the  aerodynamic  forces  acting 
on  an  aircraft  are  directly  related  to  the  wind  direction,  this  system  becomes  very  use¬ 
ful.  Specifically,  because  the  forces  imparted  on  the  aircraft  are  the  result  of  the 
aerodynamic  forces  which  are  a  function  of  the  angle  of  attack  a,  the  total  velocity 
V,  and  the  sideslip  angle  p.  For  this  axis  system,  the  origin  is  still  located  at  the  air¬ 
craft’s  center  of  gravity.  The  x-axis  is  located  in  the  direction  of  the  aircraft’  s  velocity 
vector  which  is  not  necessarily  out  the  nose.  The  y-axis  again  exits  the  aircraft’s  right 
wing,  and  the  z-axis  still  exits  the  bottom.  Due  to  the  fact  that  the  body  and  wind  sys¬ 
tem  have  the  same  origin,  the  orientation  can  be  defined  from  the  angle  of  attack  a 
and  the  sideslip  angle  p.  Components  of  the  total  velocity  V  can  be  expressed  in  terms 
of  the  body-axis  velocities  as. 


(2.1) 


u  =  VcosacosP 
V  =  VsinP 
w  =  VsinacosP 


The  vehicle  carried,  vertical-axis  system  also  has  its  origin  at  the  aircraft  center 
of  gravity.  This  axis  system  has  its  x-axis  in  the  North  direction,  the  y-axis  points  East 
and  the  z-axis  is  directed  down.  In  other  words,  the  vehicle  carried,  vertical-axis  sys¬ 
tem  is  the  earth-fixed  reference  system  only  translated  to  the  aircraft’s  center  of 
gravity.  The  method  of  comparing  this  axis  system  to  the  body-axis  system  is  by  the 
Euler  angles  \|r,  0,  and  (p  as  seen  in  Figure  2.2. 


X2,Xb  Yv 


Figure  2.2  Orientation  of  Vehicle-Carried  Vertical- Axis  to  the  Body- Axis  System 

Based  on  these  three  different  axis  systems,  the  aircraft  dynamics  can  be  de¬ 
scribed  by  12  states,  each  divided  into  4  sets  of  three  variables.  The  aircraft’s 
rotational  velocity  is  described  by  p,  q,  and  r.  The  aircraft’s  translational  motion  is 
detailed  with  V,  a,  and  p.  The  vehicle  attitude  is  found  through  \)/,  6,  and  (p.  And  the 
aircraft’s  position  is  described  by  x,  y,  and  h. 


2.1.2  Force  Equations 

With  the  different  axis  systems  described,  the  equations  of  motion  can  now  be  devel¬ 
oped  completely.  The  best  starting  place  for  this  development  is  at  the  beginning  with 
Newton’s  Second  Law. 


or 


^  =  F 
dt 


=  F 


(2.2) 


(2.3) 


In  the  above  equation,  p  is  the  linear  momentum  of  the  body,  m  is  the  total  body  mass, 
and  Vg  is  the  velocity  of  the  center  of  mass.  For  our  derivation,  mass  is  considered 
constant  since  the  fuel  weight  is  small  compared  to  the  aircraft  weight  as  a  whole. 
Thus  equation  2.3  becomes 

=  F  (2.4) 

Furthermore,  an  aircraft  is  assumed  to  behave  like  a  rigid  body,  with  angular  accel¬ 
erations.  Thus  making  it 


mi-V. 


dt  ^'OXYZ 


Ldt  ^'o'xyz 


(2.5) 


Using  the  velocity  vector  from  the  aircraft  in  the  body-axis  system,  V^.  =  ui  +  vj  +  wk 
and  the  angular  velocity  0)=  pi  +  qj  +  rk.  With  this  information,  the  force  equation 
now  become 


F^  =  m{u  +  qw-rv) 
Fy  =  miv-pw  +  ru) 
=  m{w  +  pv-qu) 


(2.6) 

(2.7) 

(2.8) 


2.1.3  Moment  Equations 

With  the  force  equations  defined  above,  the  next  step  is  to  derive  the  moment  equa¬ 
tions.  Again  the  starting  place  is  with  Newton’s  Law.  Applying  Newton’s  Laws  to  an¬ 
gular  momentum,  the  following  equation  is  developed. 


d  jj  I  _  I 

dt^OXYZ  dt^ 


o  xyz 


+  0)xH^  =  M 


(2.9) 


Given  that  Hc=Io)  with  I  being  the  inertia  matrix,  the  above  equation  becomes 


(2.10) 


M  =  /ci)  +  0)  X  /co 

Breaking  it  down  into  each  individual  axis,  the  moment  equations  are 


=  hxP-  (9-pr)  Ixy  -  (r+P^l)  hz  -  ~  hz^  (2.11) 

My  =  lyyQ  -  (/i  ”  ^ r )  "M)  "  {v^  ” ^^ )  "  4;,)/?r  (2.12) 

=  4r-  {p^-q^)hy- {q+pr)Iy^- {p-qr)I^^- (I^^-Iyy)pq  (2.13) 


2.1.4  Euler’s  Angles 

As  alluded  to  above  in  the  reference  system  section,  the  different  axis  systems  are  re¬ 
lated  through  Euler  angle  transformations.  When  transforming  coordinates  from  the 
vehicle-carried,  vertical-axis-system  to  the  body-axis  system  three  different  transfor¬ 
mations  must  occur.  The  first  is  a  rotation  about  the  Z-axis  through  an  angle  \|/.  The 
second  is  a  rotation  about  the  new  y^-axis  through  an  angle  0.  Last  is  a  rotation  about 


the  new  X2-axis  through  an  angle  <p.  At  each  rotation,  components  of  a  vector  ex¬ 
pressed  in  the  coordinate  frame  before  and  after  the  rotation  are  related  through  a  ro¬ 
tation  matrix.  Namely, 

\|r  Rotation: 


0  Rotation: 


9  Rotation: 


Xl 

cos\|/  sin\|/  0 

X 

yi 

= 

-sin\|/  cos\|t  0 

y 

h 

0  0  1 

z_ 

^2 

COS0  0  -sin0 

Xi 

yi 

= 

0  1  0 

>’1 

52 

sin0  0  COS0 

5i 

X3 

’10  o' 

X2 

>^3 

= 

0  coscp  sin  9 

^2 

^3 

0  -sin  9  cos  9 

^2 

(2.14) 


(2.15) 


(2.16) 


With  the  above  rotations  the  next  step  is  to  apply  them  in  the  specific  order  to  the  an¬ 
gular  velocity  co.  After  using  the  transformation  matrices  the  angular  velocities  can 
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be  expressed  in  terms  of  the  time  derivatives  of  the  Euler  angles, 

p  =  (p-\(tsin0  (2.17) 

q  =  0cos(p  +  \i/cos0sin9  (2.18) 

r  =  \i/cos0cos(p-0sin(p  (2.19) 

The  same  transformation  must  ^so  occur  with  the  force  on  the  aircraft  due  to  gravity. 

Applying  the  transformations  leads  to  the  following  equations. 

=  -mgsin0  (2.20) 

Fy  =  mg  cos  0  sin  9  (2.21) 

F^  =  mg  cos  0  cos  9  (2.22) 

Next,  the  linear  acceleration  terms  must  be  found.  Applying  the  transformations  will 

result  in  the  following  acceleration  equations. 

X-  ucos\(Acos0  +  v(cosv(rsin0sin<p- sinycosjp)  +  w (sinysin{p  + cos>|/sin0cos(p)  (2.23) 

y=  «sin\(fcos0  +  vcosxtfcosc(>  + vsin\(fsin0sin<p  + w  (sm>)/sm0cos(p- cosYsincp)  (2.24) 

Z  =  -«sin0  + vcos0sin(p  +  M’cos0cos((i  (2.25) 

The  final  set  of  equations  needed  are  the  angular  acceleration  equations  of  the  Euler 

angles  as  they  apply  to  the  aircraft.  As  the  aircraft  flies  its  attitudes  will  constantly  be 
changing  with  time.  This  final  set  of  equations  describe  the  state  dynamics  of  the  Eu¬ 
ler  angles  as  a  function  of  the  vehicle  angular  relationship  p,  q,  and  r. 

0  =  ^cos9-rsin9  (2.26) 

xj/  =  ^sin9sec0  +  rcos9sec0  (2.27) 

9  =  p -I- ^  sin  9  tan  0 -I- r  cos  9  tan  0  (2.28) 

This  completes  the  transformation  from  the  vehicle-carried,  vertical-axis  system  to 

the  body-axis  system  through  the  Euler  angles.  It  is  important  to  note  that  the  only 
assumptions  made  up  to  this  point  are  1)  rigid  airframe,  2)  Flat  earth  3)  Axes  fixed  to 
the  body  with  the  origin  at  the  center  of  gravity,  and  4)  Earth-fixed  reference  is  treated 
as  an  inertial  reference. 


^  2.2  Linearized  Equations  of  Motion 

Linear  aircraft  models  are  very  useful  in  aircraft  control  design.  The  quickness  in 
which  control  laws  can  be  validated  saves  many  days  of  evaluation  time  that  would 
be  spent  solely  on  a  nonlinear  simulation.  It  is  common  practice  to  first  derive  a  con- 
•  trol  law  for  the  linear  model,  and  then  apply  it  to  the  nonlinear  model.  If  the  controller 


is  robust  enough,  it  will  also  work  for  the  nonlinear  model.  Therefore,  the  next  sec¬ 
tion  will  discuss  the  derivation  of  the  linear  equations  of  motion  for  an  aircraft  from 
the  nonlinear  equations  already  presented. 

2.2.1  Linearized  Linear  Acceleration  Equations 

The  majority  of  the  forces  and  moments  applied  to  an  aircraft  are  due  to  the  aerody¬ 
namics  of  the  wings,  body,  and  tail  surfaces.  It  would  be  difficult  to  express  these 
forces  in  terms  of  u,  v,  and  w.  On  the  other  hand,  it  is  much  easier  to  express  these 
forces  and  moments  in  terms  of  vehicle  velocity  V,  angle  of  attack  a,  and  sideslip  an¬ 
gle  p.  As  shown  in  Figure  2.3,  we  can  express  the  linear  velocities  (u,v,w)  directly  in 
terms  of  V,  a,  and  P  through  the  following  relations. 

u  =  VcosPcoso  (2.29) 

V  =  Vsinp  (2.30) 

w  =  VcosPsina  (2.31) 


Figure  2.3  Linear  Velocities 


Furthermore,  equations  (2.6)-(2.8)  can  be  rewritten  as 
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u  =  rv  -  qw  -  gsinQ  + 


m 


V  =  pw  -  ru  + gsm(pcosQ  + 


w  =  qu-pv  + gcosQcos<p  + — 

m 


(2.32) 

(2.33) 

(2.34) 


The  above  linear  accelerations  can  also  be  derived  in  terms  of  V,  a,  and  P  by  differ¬ 
entiating  equations  (2-29)-(2-31)  with  respect  to  time  to  obtain 

li  =  VcosacosP- VdsinacosP- ypcosasinP  (2.35) 

V  =  VsinP  +  pVcosP  (2.36) 

V  =  VsinacosP+ VdcosacosP- V'PsinasinP  (2.37) 

Combining  the  last  two  sets  of  equations  and  expanding  them  into  separate  compo¬ 
nents,  the  following  equations  are  determined. 

V  =  (cosacosPsinG- sinPsincpcosB- sinacosPcosGcoscp)  + 


F,  F^  F- 

— cosacos3  + -2-sinp  + —  sinacosP 
m  m  m 

sinasin0  +  cosacosGcostp, 

a  =  -rsmatanp  +  ^  +  g  ( - Vcosp -  -;?cosatanp - 


(2.38) 


mVcosP 


sma  + 


cVcosP 


cos  a 


(2.39) 


P  =  ^(cosasinPsin0  +  cosPsin(pcos0-sinasinPcos6cos(p) -rcosa + 


psina- 


F  cos  a  sin  P  F  cosP  F  sin  a  sin  P 


(2.40) 


mV  mV  mV 

Now  we  begin  the  process  of  linearizing  the  above  equations.  This  is  accomplished 


by  assuming  small  perturbations  for  all  the  variables.  Specifically  the  linear  velocities 
are  now  V=Vo+AV,  a=ao-i-Aa,  and  P=Po-i-Ap.  Likewise  the  angular  velocities  are  re¬ 


placed  by  p=Po-(-Ap,  q=qo+Aq,  and  r=r(,-i-Ar.  The  force  equations  become  Fx=Fxo+ 


AFx,  Fy=Fyo+AFy,  and  Fz=F2o+AFx.  Lastly  the  airplane  attitudes  become  0=0o+A0, 
(p=(Po+A9,  and  \i/=\i/o+A\)/.  In  the  above  expression,  the  subscript  (-)o  refers  to  the 
trim  value,  and  the  A-term  is  its  perturbation.  In  order  to  finish  the  linearization,  a  spe¬ 
cific  flight  condition  must  be  chosen.  For  this  section  a  symmetric  climb  condition  is 


chosen.  This  flight  condition  has  the  following  trim  values  equal  to  zero:  Po  =  Po  = 


Qo  —  Tq  —  Fyo  —  (Po~  Vo 

Substituting  the  perturbation  equations  into  equations  (2.38)-(2.40),  the  linear¬ 
ized  equations  of  motion  can  be  developed.  The  linearization  is  done  by  neglecting 
higher-order  terms  such  as  APAa  which  is  set  equal  to  zero,  and  at  the  same  time  as¬ 
suming  small  angle  approximations  where  cosAP  =  1 ,  and  sinAP  =  Ap.  After  long 
mathematical  manipulations  the  linearized  equations  of  motion  are 

AF^cosa  AF.sina 

AV  = -^cos(0  -aJA0-i-— ^ ^  (2.41) 

Ad  =  A^-i- (-cosa^sin0^-i- sina^cos0^) - - -i- - ^ — ^2.42) 


AF  gAQ 

AP  =  A/jsina^-Arcosa^  +  ^  +  -j^cos0^ 


(2.43) 


The  next  step  is  to  convert  the  forces  from  directional  forces  in  the  Fx  and  F^  terms 
and  convert  them  to  better  known  aircraft  forces  such  as  Thrust,  Lift,  and  Drag  as  Fig¬ 
ure  2.3  shows  the  relationship  between  F^,  F^,  and  Thrust,  Lift  and  Drag. 


X  L 


Figure  2.4  Fx  and  F^  Components  in  terms  of  L,  D,  and  T 


It  is  from  this  figure  that  the  following  equations  are  derived. 


=  Lsina-Dcosa-Tcosa^  (2.44) 

F^  =  -Lcosa-Dsina- Fsina,  (2.45) 

Small  perturbation  theory  is  again  used  to  define  L=Lo+AL,  D=Do+AD,  and 

f =To+AT.  Expanding  each  term,  keeping  first-order  effects,  and  using  small  angle 
approximations,  we  obtain  a  new  set  of  equations  for  AF^  and  AF^.  These  equations 


are  then  re-substituted  into  the  linearized  equations  of  the  linear  acceleration  equa¬ 
tions  (2.41)-(2.43)  to  give  us  the  final  form  of  the  linearized  acceleration  equations. 


LAa  AD  AT 

AV  =  -ecos(6^-a  )A0-»- - +  — cos(a  -i-a) 

^  ^  o  o>  fn  mm  ^  o  t' 


(2.46) 


Ad  =  A^-|-sin(0^-a^)Ae-— Z),-  — -^sin(a^-Ha,)  (2.47) 


Aa 


AL  AT 


AB  =  Apsina  -Arcosa  +  J-cos0  Atp-i-— ^ 

t'  i'  o  o  y  o  ^  ffiY 


AY 


(2.48) 


2.2.2  Linearized  Angular  Acceleration  Equations 

With  the  linear  acceleration  terms  derived,  the  next  step  is  to  derive  the  linearized  an¬ 
gular  acceleration  equations.  Again  a  steady-level  climb  will  be  considered,  with  an¬ 
gular  velocities  redefined  as  p=Po+Ap,  q=qo+Aq,  and  r=ro+Ar.  The  moments  are  also 
defined  using  small  perturbations  about  a  trim  condition  in  order  to  make  L=Lq+AL, 
M=Mq+AM,  and  A^=Ao+AA^.  For  a  steady-level  climb,  po=qo=ro=0  and 
Lq=Mq=Nq=0.  Substituting  the  above  expressions  into  the  nonlinear  angular  acceler¬ 
ation  equations  (2.11)-(2.13)  and  using  small  angle  approximations  as  well  as  ne¬ 
glecting  higher  order  terms,  the  linearized  angular  acceleration  equations  are 

AL  =  4,Ap-4/r  (2.49) 

AM  =  lyyAq  (2.50) 

AiV  =  4,Ar-/„Ap  (2.51) 

2.2.3  Linearized  Euler  Angles 

In  order  to  linearize  the  Euler  angles,  the  perturbation  theory  applied  to  the  previous 
equations  will  be  used  again.  The  perturbations  will  be  applied  to  equations  (2. 17)- 


(2.19)  with  the  same  simplified  assumptions  of  small  angle  approximation  and  dele¬ 
tion  of  higher  order  terms.  Making  these  substitutions  and  applying  the  small  pertur¬ 
bation  assumptions  yields  the  following  equations 


A0  = 

(2.52) 

A  •  _ 

COS0 

0 

(2.53) 

A9  =  A/?-l-  tan0^Ar 

(2.54) 

2.3  Nonlinear  Simulation  Model 

The  twelve  nonlinear  simultaneous  second-order  differential  equations  derived  in  the 
previous  sections  are  implemented  in  the  simulation  in  the  form  of  a  MATLAB  spe¬ 
cial  function,  or  s-function.  The  exact  methodology  of  the  MATLAB  function  is  not 
important.  What  is  important  is  the  basic  format  of  the  function,  as  well  as  its  inte¬ 
gration  with  the  simulation  program  SIMULINK.  The  next  sections  will  outline  these 
points. 

2,3.1  Component  Integration-4he  S-Function 

Systems  that  are  saved  in  the  SIMULINK  window  as  s-functions  are  treated  in  a  sim¬ 
ilar  manner  as  MATLAB  treats  m-functions.  In  MATLAB,  the  program  calls  a  certain 
subprogram  to  be  run.  Likewise,  SIMULINK  calls  a  certain  subprogram  to  be  run. 
The  subprogram  must  be  in  the  same  directory,  but  that  is  the  only  stipulation  on  run¬ 
ning  it.  The  s-function  makes  it  possible  for  the  user  to  write  specific  routines  for  the 
simulation.  The  F-15  dynamics  that  were  acquired  from  Wright-Patterson  AFB  are 
located  in  two  primary  routines,  f25aero  and  f25eng,  quantifying  the  aircraft  aerody¬ 
namic  and  propulsion  characteristics  respectively.  The  routine  atmos  calculates  the 
current  atmospheric  parameters  such  as  density  and  temperature  to  be  used  by  the 
previous  two  routines.  All  three  of  these  routines  were  provided  in  FORTRAN  format 
from  the  AIAA  Design  Challenge  [5].  The  s-function  used  for  the  open-loop  simula¬ 
tion  is  shown  in  Appendix  A. 

Initial  conditions  for  the  aircraft  states  are  also  established  inside  the  s-function. 
These  particular  values  are  very  important  when  determining  the  aircraft  trim  points 
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for  a  given  flight  condition.  The  s-function  also  makes  it  possible  to  determine  the 
output  variables.  For  the  open-loop  case  the  12  states  were  chosen  in  an  effort  to  val¬ 
idate  the  equilibrium  state  responses  at  the  two  chosen  flight  conditions. 

2.3.2  The  SIMUIjJNK  Model 

As  far  as  the  SMULINK  model  is  concerned,  its  job  is  made  easier  by  the  s-function. 
SIMULINK  only  needs  to  set  up  the  inputs  and  outputs  of  the  s-function.  A  clock  is 
also  needed  in  order  to  have  a  time  response  for  the  actual  simulation.  That  is  all  there 
is  to  the  SIMULINK  model.  As  seen  in  Figure  2.5  the  inputs  to  the  s-function  are  cho¬ 
sen  based  on  the  inputs  needed  to  run  the  s-function.  Likewise  the  outputs  were  de¬ 
termined  based  on  the  outputs  deemed  necessary  by  the  user  when  he  creates  the  s- 
function. 
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Figure  2.5  SIMULINK  model 


The  inputs  for  the  s-function  f25sfn  are  typical  inputs  into  most  aircraft.  These  in¬ 
clude  the  horizontal  stabilizers  both  differential  (DD)  and  symmetric  (DH).  Ailerons 
(DA)  and  rudder  (DR)  are  also  needed.  And  finally,  the  thrust  from  the  engines  is  also 
an  input.  The  inputs  PLAR  and  PLAL  refer  to  the  right  and  left  power  level  angles 
respectively,  PLASYM  assures  that  both  are  set  equally.  As  for  the  outputs,  they  were 
chosen  to  be  the  12  states. 

2.4  Linearization  of  the  Model 

It  has  already  been  discussed  that  a  linear  model  facilitates  the  design  process.  Grant¬ 
ed  several  different  flight  points  must  be  designed  for  in  order  to  use  the  controller  in 
the  nonlinear  environment.  Therefore,  different  flight  points  must  be  chosen.  In  pre¬ 
vious  work,  Lt.  James  Dutton  chose  two  flight  points  in  order  to  set  up  his  control  law 
[10].  In  order  to  further  his  research  and  to  validate  my  own,  I  chose  to  use  the  same 
flight  points.  The  original  decision  for  the  flight  points  was  not  just  a  guess.  They 
were  chosen  through  the  use  of  a  program  that  calculates  the  equilibrium  states  of  an 
aircraft,  namely  trim.  This  program  doesn’t  allow  any  variable  to  remain  constant. 
Instead  it  varies  the  trim  values  in  an  effort  to  find  an  equilibrium  point.  By  allowing 
these  variations,  the  two  trim  points  were  determined  for  this  report.  Table  2.1  shows 
how  the  initial  guess  settled  at  the  trim  point  for  flight  point  #1. 


Table  2.1:  TWin  Data  FPl  (9,800ft,  0.5  M) 


Parameter 

Units 

Initial 

Guess 

Trim  Value 

State 

Derivative 

V 

ft/s 

539.08 

539.08 

r7.1636-^ 

a 

rad 

0.0801 

0.0801 

-1.8594-0 

q 

rad/s 

0 

.2,79-19 

-2.1328-^ 

0 

rad 

0.0801 

0.0803 

-2.7949-10 

p 

rad/s 

0 

-2.54-20 

2.3895-10 

<p 

rad 

0 

-1.02-^* 

-2.8053-20 

r 

rad/s 

0 

-3  22-20 

-1.8936-20 

¥ 

rad 

0 

-1.66-^^ 

-3.2388-20 

P 

rad 

0 

-8.41-21 

-2.9103-20 

5h 

deg 

-2.827 

-2.880 

N/A 

5pLA 

deg 

37.4 

37.4 

N/A 

5a 

deg 

0 

0 

N/A 

5r 

deg 

0 

0 

N/A 

Y 

rad 

0 

0.0002 

N/A 

Once  the  two  flight  points  have  been  determined,  linear  model  describing  the  F- 15  at 
a  specific  flight  point  can  be  obtained  in  the  form  of  the  following  matrix  equations. 

x(t)  =  Ax(t)  +Bu(t)  (2.55) 

y(0  =  Cx(t) +Du(t)  (2.56) 

The  data  collected  from  the  trim  program  was  converted  into  the  system  of  equations 

above,  and  then  saved  under  two  different  files  to  be  used  later.  The  first  flight  point 
is  at  9,800  ft  flying  at  Mach  0.5  which  translates  to  539.0  ft/sec.  The  second  flight 
point  is  at  30,000  ft  again  flying  at  Mach  0.5.  For  this  altitude,  the  velocity  is  497  ft/ 
sec.  The  files  that  contain  the  matrix  information  are  trim539ss  and  trim497ss  for 
flight  points  FPl  and  FP2  respectively.  These  two  files  were  created  by  the  MATLAB 
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function  linmod,  which  took  the  aircraft  data  and  translated  it  into  the  respective  state 
space  models. 


# 


#  Chapter  3  F-15  Open-Loop  Analysis 

The  purpose  of  this  chapter  is  to  discuss  the  open-loop  analysis  of  the  F-15.  Open 
loop  analysis  is  very  important  at  the  initial  stage  of  a  control-law  design.  If  an  air- 

0  craft  is  inherently  unstable,  it  will  manifest  itself  in  the  open-loop  simulation.  This 

will  in  turn  provide  engineers  with  knowledge  concerning  which  modes  will  need 
what  type  of  modifications  in  order  to  control  the  aircraft. 

Open-loop  analysis  is  also  important  to  linear  simulation  verification.  It  has  been 

•  previously  stated  that  linear  design  points  are  used  by  engineers  in  order  to  design  an 

adequate  control  law.  If  responses  from  the  linear  simulation  and  the  nonlinear  sim¬ 
ulation  are  very  different,  the  control  law  designed  for  the  nonlinear  simulation 

^  simply  will  not  work.  Therefore,  by  looking  at  the  open-loop  time  responses  of  the 

F-15,  the  designer  will  know  whether  or  not  the  linear  simulation  accurately  de¬ 
scribes  the  nonlinear  case. 


3.1  Linearized  Model 

The  nonlinear  model  for  the  F-15  has  already  been  shown  in  detail.  In  order  to  run  a 
similar  simulation  for  the  linear  case  one  must  create  a  separate  SIMULINK  file.  In 
Figure  3.1  one  sees  the  SIMULINK  file  used  to  test  the  linear  case. 
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Notice  that  the  inputs  and  outputs  are  exactly  the  same  as  in  the  nonlinear  case  de¬ 
picted  in  Figure  2.5  for  consistency.  As  long  as  the  inputs  are  the  same  and  the  outputs 
measured  are  the  same,  it  stands  to  reason  that  the  linear  model  and  the  nonlinear 
model  have  a  good  chance  of  having  the  same  time  responses. 

3.2  Evaluation  of  the  Linearized  Model 

The  nonlinear  and  linear  models  were  run  with  a  variety  of  inputs  in  order  to  validate 
the  two  systems.  As  you  will  see  in  the  next  few  pages,  the  linear  simulation  accu¬ 
rately  describes  the  nonlinear  dynamics  found  in  the  F-15.  For  most  cases,  nonlinear 
model  had  more  fluctuations,  but  the  overall  responses  of  the  linear  simulations  were 
accurate  enough  to  be  used  for  the  control-law  application. 

3.2.1  Short-Period  Time  Response 

The  short  period  mode  was  excited  by  inputting  a  20-second  pulse  of  2  degrees 
down  on  the  elevator  of  the  aircraft  at  both  flight  points.The  aircraft  time  responses 
are  found  in  Figure  3.2  through  Figure  3.5.  In  these  figures,  the  dotted  line  represents 
the  linear  response  and  the  solid  line  is  the  nonlinear  time  response. 

3.2.2  Roll  Time  Response 

For  Flight  Point  1,  a  pulse  was  applied  to  the  right  aileron  down  one  degree  for 
20  seconds  in  order  to  check  the  roll  mode  of  the  aircraft.  Again  the  linear  simulation 
closely  matched  the  nonlinear  case.  For  Flight  Point  2,  the  right  aileron  was  deflected 
down  2  degrees  for  20  seconds.  The  linear  and  nonlinear  time  responses  are  compa¬ 
rable.  These  time  responses  can  be  found  in  Figures  3.6-3.9,  with  the  dotted  lines 
again  being  the  linear  case  and  the  solid  lines  being  the  nonlinear  case. 

3.2.3  Dutch-Roll  Time  Response 

For  both  flight  points,  the  rudder  was  deflected  1  degree  to  the  left  for  20  seconds 
in  order  to  check  the  Dutch  roll  mode  of  the  F-15.  Both  the  nonlinear  case  and  the 
linear  case  had  similar  time  responses.  In  Figures  3.10-3.13,  the  time  responses  are 
shown  with  the  dotted  line  depicting  the  linear  case  and  the  solid  showing  the  non- 


linear  case. 
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Figure  3.3  Aircraft  Responses  to  a  20-second  Elevator  Pulse  of  2°. 
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Figure  3.5  Aircraft  Responses  to  a  20-secon(i  Elevator  Pulse  of  2°. 


Figure  3.6  Aircraft  Responses  to  a  20-second  Aileron  Pulse  of  1®. 
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Figure  3.7  Aircraft  Response  to  a  20-second  Aileron  Pulse  of  1°. 


Flight  Point  #2  Open  Loop  Aileron  Pulse  2  deg 


Figure  3.8  Aircraft  Response  to  a  20-second  Aileron  Pulse  of  2° 
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Flight  Point  #2  Open  Loop  Aileron  Pulse  2  deg 


time  (sec) 


Figure  3.9  Aircraft  Response  to  a  20-second  Aileron  Pulse  of  2°. 


Flight  Point  #1  Open  Loop  Rudder  Pulse  1  deg 


Figure  3.10  Aircraft  Response  to  a  20-second  Rudder  Pulse  of  1°. 


Figure  3.1 1  Aircraft  Response  to  a  20-second  Rudder  Pulse  of  1°. 


Flight  Point  #2  Open  Loop  Rudder  Pulse  1  deg 


Figure  3.13  Aircraft  Response  to  a  20-second  Rudder  Pulse  of  1°. 


Chapter  4  Lateral  Control  Using  TECS 

4.1  Background. 

Modem  control  theory  have  focused  on  the  improvement  of  the  controller  perfor¬ 
mance  and  robustness  to  plant  uncertainties.  Such  multivariable  control  synthesis  in¬ 
cludes  linear  quadratic  Gaussian/loop  transfer  recovery  (LQG/LTR),  H2,  and 
synthesis.  However,  these  techniques  are  of  high  order,  and  hence  not  very  cost  ef¬ 
fective.  In  order  to  lower  the  cost,  classical  design  procedures  based  on  single-loop 
closure  and  root  locus  have  been  applied  to  multi-loop  control  systems.  However, 
these  designs  can  have  low  performance  and  poor  robusmess  to  model  uncertainties. 
The  trade-off  between  performance  and  robustness  is  not  a  simple  process  and  can  be 
very  time  consuming.  Furthermore,  single-loOp  analysis  does  not  incorporate  the 
multivariable  aspects  of  a  design.  This  fact  will  also  limit  the  maximum  achievable 
performance.  Single-loop  design  is  easier,  but  not  as  effective  nor  as  robust.  Multi- 
variable  design  would  have  better  performance  and  robustness,  but  is  of  high  order. 

The  TECS  controller  is  a  combination  of  the  simplicity  of  the  single-loop  design, 
linked  with  the  performance  and  robustness  of  the  multivariable  methods.  Single¬ 
loop  control  analysis  falls  short  where  cross-coupling  effects  are  concerned.  For  in¬ 
stance,  most  single-loop  controller  have  a  separate  loop  for  the  throttle  and  the 
elevator.  This  results  in  either  over  or  under  control  of  the  surfaces  in  order  to  per¬ 
form  flight  maneuvers.  The  TECS  controller  works  by  converting  energy  states.  It 
positions  the  control  surfaces  simultaneously  in  order  to  achieve  the  maximum  per¬ 
formance  with  little  or  no  conflicting  commands  from  the  control  surfaces.  The  exact 
manner  is  described  in  the  next  section. 

4.2  Development  of  the  TECS  Concept 

As  stated  above,  the  TECS  controller  functions  by  altering  the  energy  states  of  the 
aircraft.  Specifically,  this  means  trading  off  potential  or  kinetic  energy  in  order  to  per¬ 
form  a  specific  flight  function.  The  key  equations  governing  the  total  energy  control 
concept  were  developed  in  [2  4  7  8].  They  are  reviewed  here  for  completeness  and  to 
define  the  structure  of  the  control  system. 
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The  total  energy  E(t)  of  the  aircraft  treated  as  a  point  mass  is  given  by 

E  it)  =  (t)  (t)  gh  it)  (4.1) 

where  V(t)  is  the  aircraft’s  total  velocity,  m(t)  is  the  mass  of  the  aircraft,  h(t)  is  the 
altitude,  and  g  is  the  gravitational  acceleration.  The  total  energy  rate  is  then  given  by 
differentiating  equation  (4.1)  with  respect  to  time  and  assuming  small  flight  path  an¬ 
gle  and  that  the  aircraft  mass  is  constant  or  slowly  varying. 

E'(t)  =  mit)gVit)  i'^  +  yit))  (4.2) 

Now,  the  equation  can  be  rewritten  in  terms  of  thrust  required,  T^gq.  Using  the  equa¬ 
tions  of  motion  developed  in  Chapter  2  we  obtain 

T^^git)  =  mit)gi^+yit))  +Dit)  (4.3) 

o 

Assuming  that  the  thrust  is  trimmed  against  the  drag,  and  that  the  variation  in  drag  is 
slow,  then  the  change  in  thrust  required  is  proportional  to  the  change  in  energy  divid¬ 
ed  by  the  aircraft  velocity.  This  in  turn  leads  directly  to 


=  m(r)g(-^-hY(0)  (4.4) 

o 

This  last  equation  is  important  because  it  shows  that  the  total  energy  of  the  aircraft 
can  be  regulated  directly  using  thrust.  In  essence,  the  energy  level  of  the  aircraft  is 
controlled  by  the  throttle.  It  also  shows  that  the  distribution  of  the  energy  cannot  be 
controlled  by  the  throttle,  hence  the  need  for  the  elevator.  The  elevator  distributes  the 
energy  available  between  kinetic  and  potential  energies. 

With  the  two  control  surfaces  having  their  defined  roles,  integral  and  proportional 
control  can  be  applied  to  the  throttle  and  elevator.  The  following  equations  result 
from  this  application  with  5t  and  6e  being  the  throttle  and  elevator  commands  respec¬ 


tively. 


KtI  (V: 


8,  =  mg(K^p+-^) 

u 


\8 


+  Yp 


5,  =  + 


K 


El 


)\j-\ 


(4.5) 


(4.6) 


where  Ye=Y'Yc  andVe  =  V  -  represent  the  errors  of  the  flight  path  angle.  Ye.  and 


acceleration,  Vg  from  the  commanded  values  and  respectively.  Implementation 
of  this  system  into  SIMULINK  format  is  shown  in  Figure  4.1. 


Figure  4.1  TECS  Controller  Structure  for  the  Longitudinal  Axis 
For  the  above  block  diagram,  the  parameters  K-pp  K-n,  Kep^  and  Kej  are  the  propor¬ 
tional  and  integral  feedback  gains  on  the  throttle  and  the  elevator  respectively.  The 
gains  Kq  and  Kq  are  likened  to  a  stability  augmentation  system  (SAS)  that  improves 
the  short  period  response.  The  gains  Ky,  and  Kjj  provide  feedback  correction  to  the 
airspeed  and  altitude  errors.  Lastly,  the  gains  and  are  scheduled  accord- 


ingly  to  compensate  for  the  calibrated  airspeed  and  the  fluctuations  in  gross  weight 
respectively.  This  completes  a  cursory  explanation  of  the  TECS  controller  structure 
for  the  longitudinal  axis.  This  report  is  concerned  with  the  lateral  dynamics  of  the 
McDonnell  Douglas  F-15,  flierefore  the  next  section  will  detail  the  TECS  controller 
structure  for  the  lateral  axis.  .  ^ 

4.3  Lateral  TECS  Controller 

The  total  energy  conservation  principle  is  not  as  clearly  seen  in  the  lateral  axis  as  it 
is  in  the  longitudinal  axis.  For  the  lateral  axis  the  energy  balance  concerns  a  coordi¬ 
nated  turn.  As  an  aircraft  starts  to  turn,  the  natural  tendency  is  for  the  nose  to  yaw 
against  the  direction  of  the  turn,  causing  a  slip  through  the  air.  This  is  not  optimum 
performance  of  an  aircraft.  In  order  to  perform  the  turn  in  a  coordinated  manner,  the 
pilot  needs  to  add  rudder  in  order  to  counter  the  adverse  yaw  due  to  aileron  deflection 
and  the  induced  drag  causing  the  yaw.  The  lateral  TECS  controller  has  as  its  founda¬ 
tion  an  energy  balance  between  rolling  motion  (rolling  energy)  and  yawing  motion 
(yawing  energy)  in  order  to  perform  a  coordinated  turn.  The  side  force  equation  and 
the  requirement  for  a  turn-coordination  is  depicted  in  equation  (4.7). 

pw  (4.7) 

For  a  coordinated  turn  in  steady  state,  the  aircraft  bank  angle  and  the  yaw  rate  must 
have  a  certain  relationship.  This  relationship  is  shown  in  Equation  4.8. 

ts-Y- 

Breaking  the  maneuver  into  small  incremental  steps  is  very  helpful.  The  change  in 
bank  angle  is  then  given  by 

■ 

A(t)  =  —  (p-hAr)  (4.9) 

s 

Motion  in  the  longitudinal  axis  movement  is  mainly  governed  by  the  throttle  and  the 
elevator.  The  lateral  axis  dynamics  are  controlled  by  the  aileron  and  the  rudder. 
Therefore,  it  only  makes  sense  to  have  an  aileron  command  that  is  a  function  of  Atp, 

and  have  the  rudder  be  a  function  of  (P  -  Ar) .  When  this  is  done,  the  rudder  func- 
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tions  effectively  as  a  weathervane,  so  that  at  a  constant  bank  angle,  the  difference 
(P  -  Ar)  is  equalled  to  zero. 

Proportional  and  integral  control  are  used  for  consistency  in  order  to  develop  ad¬ 
equate  control  responses  for  the  TECS  controller.  Applying  proportional  and  integral 

control  on  the  quantities  ((3  -  Ar)  and  (3  +  Ar) ,  we  develop  the  aileron  and  rudder 
command  in  equations  4.10  and  4. 1 1 . 

5„  =  Kr(Kyf+  (4.11) 

IJ 

Implementation  of  the  previous  equations  into  SIMULENK  is  found  in  Figure  4.2. 


Figure  4.2  TECS  Controller  Structure  for  the  Lateral  Axis 
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The  gains  in  lateral  TECS  controller  are  directly  comparable  to  the  longitudinal  gains 
described  in  the  previous  section.  Specifically  the  parameters  K^p,  Kyp  Krj  and  Kyj 
are  the  proportional  and  integral  gains  on  the  aileron  and  rudder  surfaces.  The  gains 
Kp  and  Kq,  are  used  for  the  roll  control  system  while  is  incorporated  in  a  yaw 

damper  to  better  the  Dutch-roll  characteristics.  The  parameters  BCp  and  Ky  are  pro¬ 
portional  gains  on  the  sideslip  and  heading  errors,  and  finally  and  Kj.  are  used  for 
gain  scheduling  for  the  appropriate  flight  dynamic  pressures. 

The  above  gains  will  be  determined  by  a  design  method  based  on  parameter  op¬ 
timization,  namely  the  SANDY  program.  This  optimization  process  requires  a  cost 
function  that  takes  into  account  the  command  tracking  problem  as  well  as  disturbance 
attenuation  and  good  closed-loop  stability.  One  such  cost  function  is  defined  as, 

J=  Urn  f£[0i(p(r)-P^(O)^  +  (22(V(O-Ve(O)^]  (4.12) 

tf 

where  Pc(t)  =  Pcmd(l"®'^^)  Vc(0=¥cmd(l'®'^^)-  parameters  Qj  and  Q2  are  se¬ 

lected  to  apply  the  proper  weights  to  create  trade-offs  on  the  tracking  variables.  An¬ 
other  added  constraint  to  this  cost  function  is  to  establish  guidelines  for  the 
optimization  process.  Specifically  a  certain  eigen  value  and  damping  ratio  were  cho¬ 
sen  as  follows, 

RealiK.)  =  -0.1  (4.13) 

C,>L;„  =  0.7  (4.14) 

The  m-file  zzlatdes.m  was  used  to  develop  the  necessary  data  and  put  it  in  the  proper 

format  for  the  program  SANDY.  A  detailed  discussion  of  this  process  is  located  in 
the  next  chapter. 


# 


9 


Chapter  5  TECS  Linear  and  Nonlinear  Evaluations 

The  focus  of  this  chapter  is  to  first  discuss  the  design  process  for  the  TECS  con¬ 
trol  law  starting  with  the  linear  model,  then  report  the  results  on  the  linear  model. 
Next,  the  gains  acquired  for  the  linear  model  will  be  tested  on  the  nonlinear  simula¬ 
tion.  After  testing  the  gains  I  will  provide  an  evaluation  of  the  overall  design  process 
as  well  as  the  results. 

5.1  Linear  Closed-Loop  Model  Design  and  Evaluation 

The  open-loop  evaluation  was  covered  in  Chapter  3,  therefore  this  section  will 
only  concern  itself  with  the  closed-loop  system. 

5.1.1  Initial  Model  Synthesis 

Several  tools  were  available  in  order  to  create  an  optimum  solution  to  the  closed-loop 
linear  model  with  a  TECS  control  law  for  the  lateral  dynamics  of  the  McDonnell 
Douglas  F-15.  The  major  question  was  whether  or  not  I  would  use  the  tools  correctly. 
The  primary  tool  was  a  numerical  optimization  program,  SANDY.  This  program  with 
a  given  set  of  inputs  would  numerically  calculate  the  optimum  gains  for  the  control 
law.  Before  the  gains  could  be  calculated  the  system  had  to  be  properly  input  into  the 
SANDY  program.  It  only  stands  to  reason  that  SANDY  like  most  computer  software 
only  works  as  well  as  you  set  it  up  to  work.  The  old  adage  “garbage  in  garbage  out” 
is  very  true  to  the  working  of  SANDY.  In  previous  course  applications,  I  was  given 
a  file,  latdes.m,  that  would  send  SANDY  the  proper  information  in  the  proper  format 
for  the  lateral  dynamics  of  a  TECS  control  law.  The  only  difference  was  that  this  file 
dealt  in  degrees  and  per  g  units  and  it  was  designed  for  the  Boeing  767  and  not  a  F- 

15.1  used  latdes.m  as  a  template  in  order  to  create  zz]atdes.m  which  put  linear  model 
of  the  F-15  in  the  correct  units  and  in  the  proper  order  for  SANDY  to  use.  The  output 
of  this  file  was  zzlatTECSl.dat  and  zzlatTECS2.dat  corresponding  to  flight  point  1 
and  flight  point  2  respectively.  These  MATLAB  files  can  be  found  in  Appendix  C. 
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5.1.2  SANDY  Design 

The  two  above  files  were  implemented  into  the  SANDY  program.  Initially,  SANDY 
was  set  only  to  run  for  a  short  time.  Like  most  numerical  optimization  programs,  this 
can  be  a  lengthy  process.  By  limiting  the  number  of  iterations,  problems  could  be 
found  and  corrected  before  they  compounded  into  an  unfavorable  solution.  It  was  a 
good  thing  the  number  was  limited  because  several  unit  miscalculations  were  found 
and  fixed  before  too  much  computation  time  accrued.  Once  all  the  corrections  were 
made,  the  gains  for  the  second  flight  point  came  very  quickly  with  acceptable  results. 
Design  of  flight  point  1  was  a  little  more  difficult  to  achieve.  Specifically,  the  gain 
and  phase  margins  were  unacceptable  for  the  aileron  command  loop.  The  gain  margin 
was  -8.894  dB  and  the  phase  margin  was  -29.21  degrees.  The  design  created  for  the 
linear  system  will  have  to  be  implemented  on  a  nonlinear  system,  thus  robustness  of 
the  controller  needed  improvement.  As  a  guide  I  used  ±6  dB  and  ±45°  for  the  gain 
and  phase  margins  respectively. 

In  order  to  correct  this  shortcoming,  the  zzlatTECSl.dat  file  was  modified.  The 
aileron  was  using  too  much  control,  thus  to  improve  the  robustness  I  limited  the  aile¬ 
ron  authority  by  placing  an  upper  limit  on  its  use.  I  needed  to  lower  the  upper  limit 
several  times  in  order  to  get  acceptable  results.  These  results  are  located  in  Table  5.1 
which  shows  the  worst-case  eigenvalues  and  Table  5.2  which  has  the  gain  and  phase 
margins  for  both  controllers  at  both  flight  points.  The  gains  calculated  by  SANDY 
are  found  in  Appendix  D. 


Table  5.1:  Closed-Loop  Stability  Characteristics 


Flight  Point 

Eigenvalues 

Damping 

Frequency 

1 

-0.1 

1.0 

0.1 

-  4.49  ±  4.54i 

0.7 

6.357 

2 

-0.1 

1.0 

-  0.326  ±  0.333/ 

0.7 

Table  5.2:  Single-Loop  Stability  Margins 


Flight  Point 

5aG.M. 

5a  P.M. 

5rG.M. 

5rP.M. 

1 

-37.13 

-47.39 

OO 

-65.14 

2 

OO 

-82.27 

OO 

-70.72 

These  results  were  calculated  by  a  file  titled  sandybl.m.  This  file  would  take  the 
output  from  the  SANDY  program  and  evaluate  each  command  loop  of  the  closed- 
loop  system  by  breaking  the  respective  control  loop.  The  Bode  plots  of  the  closed- 
loop  system  for  each  flight  point  and  each  control  surface  are  in  Figure  5.1  though 
Figure  5.4. 


Frequency  (Rad/Sec) 


Figure  5.1  Flight  Point  1  Aileron  Bode  Plot 
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Rudder  Bode  Plot  For  Flight  Point  2 


5.1.3  Linear  Model  Closed-Loop  Command  Responses 

The  closed-loop  command  responses  were  determined  through  the  use  of  the  m-file 
linan.m  which  takes  the  aircraft’s  closed-loop  matrices  and  inputs  them  into  a  linear 
simulation  function  of  MATLAB.  The  response  for  the  sideslip  p,  heading  \)/,  aileron 
5a,  and  rudder  5r  are  discussed  below.  In  all  figures  describing  the  time  responses  the 
thicker  line  is  the  input  command  and  the  thin  line  is  the  MATLAB  linear  solution. 
Heading  Command 

The  linearized  aircraft  responses  to  a  30-degree  heading  command  are  shown  for 
both  flight  conditions  in  Figures  5.5  and  5.6  respectively.  The  heading  tracking  is 
very  good  in  both  cases  showing  no  overshoot  and  very  few  deviations  from  the  com¬ 
mand.  In  addition,  the  amount  of  rudder  and  aileron  deflections  to  acquire  these  re¬ 
sults  are  well  within  component  limits  as  well  as  within  common  sense  parameters. 
Sideslip  Command 

The  aircraft  was  subjected  to  a  5-degree  constant  sideslip  command  in  order  to  eval¬ 
uate  the  sideslip  performance.  Flight  point  2  tracked  the  sideslip  command  very  well 


with  no  overshoots.  In  addition,  the  control  deflections  are  moderate  without  large 
fluctuations  that  would  add  vibrations  to  the  aircraft  For  flight  point  1,  the  aircraft 
did  not  track  closely  the  command.  However,  the  time  response  was  reasonable 
enough  to  continue  the  design  process.  The  control  deflections  for  this  command  are 
also  within  tolerances.  Figures  5.7  and  5.8  show  the  sideslip  command  responses  for 
flight  points  1  and  2  respectively. 


time 


Figure  5.5  Flight  Point  1  Linear  Aircraft  Response  to  a  30°-Heading  Command 


Flight  Point  2  Heading  Command  =  30  degs 


Figure  5.6  Flight  Point  2  Linear  Aircraft  Response  to  a  30 -Heading  Command 


Flight  Point  1  Sideslip  Command  =  5  degs 


Figure  5.7  Flight  Point  1  Linear  Aircraft  Response  to  a  5 -Sideslip  Command 
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Flight  Point  2  Sideslip  Command  =  5  degs 


time 

Figure  5.8  Flight  Point  2  Linear  Aircraft  Response  to  a  5°-Sideslip  Command 
5.1.4  Evaluating  the  Linear  SIMULINK  Model 

Up  until  this  point  the  linear  simulation  of  the  F- 15  was  done  in  the  MATLAB  envi¬ 
ronment.  The  nonlinear  simulation  is  set  up  in  the  SIMULINK  environment.  Previ¬ 
ous  assignments  have  proven  that  what  is  in  the  MATLAB  environment  will  not 
always  directly  transfer  into  the  SIMULINK  environment  correctly.  Therefore,  the 
linear  SIMULINK  environment  must  be  tested  and  the  results  verified  with  those  at¬ 
tain  from  the  MATLAB  functions. 

The  first  step  of  the  process  was  to  extract  the  gains  acquired  from  the  SANDY 
program.  The  programs  zzlatTECSlc.m  and  zzlatTECS2c.in  were  both  modified 
with  the  gains  explicitly  stated  so  the  SIMULINK  program  would  read  them  correct¬ 
ly.  The  next  modification  came  from  a  units  conversion.  The  states  in  the  aircraft 
model  are  all  in  radians.  SANDY  solved  for  the  states  in  degrees,  and  the  inputs  into 
the  aircraft  model  are  also  in  degrees.  In  order  to  keep  the  units  consistent,  the  outputs 
of  the  state-space  model  needed  to  be  converted  to  degrees.  This  was  done  in  the 
trim497ss.m  and  trim539ss.m  files  before  they  were  entered  into  the  SIMULINK 


representation.  These  two  MATLAB  files  can  be  found  in  Appendix  B. 

The  SIMULINK  model  consists  of  both  the  lateral  and  the  longitudinal  models, 
therefore  the  gains  gotten  by  Lt.  Dutton  needed  to  be  added  before  the  simulation 
could  be  run.  Lt.  Dutton  saved  his  gains  under  f25sd539_cx.inat  and  f25sd497_cx.- 
mat  for  flight  points  1  and  2  respectively  [10].  The  complete  SIMULINK  model  is 
located  in  Figure  5.9. 


Figure  5.9  Linear  SIMULINK  Model  Representation 


Due  to  the  fact  that  this  was  a  linear  model,  the  numerical  integration  method 
used  to  propagate  the  simulation  was  linsim  which  is  found  within  the  simulation  pa¬ 
rameters  menu.  A  few  bugs  were  discover  during  this  process  and  they  were  quickly 
corrected.  In  the  end,  the  results  of  the  linear  simulation  provided  by  MAT  ALB  were 
consistent  with  the  simulation  provided  by  the  linear  simulation  done  in  the  SIM- 
ULINK  environment  The  results  of  this  comparison  are  found  in  Figures  5.10-5.13. 
The  heavier  line  is  the  combination  of  the  MATLAB  and  SIMULINK  solutions  while 
the  thin  line  is  the  original  command  input  that  the  simulations  are  trying  to  track. 


SIMULINK  vs  Mallab  Comparison  FP#1 


Psicmd  =  30  degs 


time  (sec) 


Figure  5.10  Flight  Point  1  Linear  Comparison  Heading  Command  30®. 


Betacmd  =  5  degs 


SIMULINK  vs  Matlab  Comparison  FP#1 


Figure  5.1 1  Flight  Point  1  Linear  Comparison  Sideslip  Command  5°. 


Figure  5.12  Flight  Point  2  Linear  Comparison  Heading  Command  30°. 
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SIMULINK  vs  Matlab  Comparison  FP#2  Betacmd  =  5  dags 


time  (sec)  time  (sec) 

Figure  5.13  Flight  Point  2  Linear  Comparison  Sideslip  Command  5°. 

5.2  Nonlinear  Closed-Loop  Model  Evaluations 

In  order  to  do  the  nonlinear  closed-loop  evaluations,  the  gains  acquired  from  the  lin¬ 
ear  closed-loop  calculations  were  inserted  into  the  nonlinear  simulation  shown  in 
Chapter  2  of  this  report.  The  command  responses  were  repeated  for  the  nonlinear 
simulation  with  favorable  results. 

5.2.1  Nonlinear  Closed-Loop  Command  Responses 

The  closed-loop  time  responses  were  calculated  through  the  use  of  the  SIMULINK 
model.  The  SIMULINK  model  has  command  inputs  to  be  defined  by  the  user.  The 
previous  commands  were  a  5-degree  sideslip  command  and  a  30-degree  heading 
change,  each  done  separately.  These  commands  were  repeated  for  the  nonlinear  mod¬ 
el. 
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Heading  Command 

The  nonlinear  simulation  was  tasked  to  track  a  30-degree  heading  change  for  both 
flight  points.  For  flight  points  1,  the  nonlinear  and  linear  cases  produced  nearly  iden¬ 
tical  results.  The  nonlinear  solution  tracked  the  heading  command  very  accurately 
with  no  overshoots.  In  addition,  the  control  deflections  were  well  within  tolerable 
limits.  One  difference  between  the  linear  and  nonlinear  control  deflections  was  there 
was  a  single  fluctuation  in  control  surfaces  found  in  the  nonlinear  simulation  that  was 
not  present  in  the  linear  case.  This  only  makes  sense  that  noise  in  the  control  surfaces 
would  represent  itself  in  the  nonlinear  case  and  not  the  linear  case.  Flight  point  2  also 
tracked  the  heading  command  very  well.  As  with  the  first  flight  point,  flight  point  2 
had  no  overshoots  and  quickly  matched  the  input  command.  The  control  deflections 
for  this  flight  point  had  more  fluctuations  when  compared  to  the  linear  case.  However, 
the  deflections  consistently  matched  the  linear  case  deflections.  Results  of  the  simu¬ 
lations  were  graphed  simultaneously  in  Figures  5. 14-5. 17  in  order  to  show  the  differ¬ 
ences  described  above  graphically.  In  general,  the  thin  line  represents  the  linear 
model  whereas  the  thicker  line  represents  the  nonlinear  model.  The  exception  is  on 
the  first  graph  of  each  set,  where  the  input  command  is  also  plotted.  In  these  graphs 
the  dashed  line  is  the  linear  case,  the  thick  line  is  the  nonlinear  case,  and  the  thin  line 
is  the  command  input. 


Psicmd  =  30  degs 


Nonlinear  vs  Linear  Time  Response  tor  FP#1 


Figure  5.14  Flight  Point  1  Nonlinear  vs.  Linear  Time  Responses 


Nonlinear  vs  Linear  Time  Response  tor  FP#1 


Psicmd  =  30  degs 


Nonlinear  vs  Linear  Time  Response  for  FP#2 


Figure  5.16  Flight  Point  2  Nonlinear  vs.  Linear  Time  Responses 


Nonlinear  vs  Linear  Time  Response  for  FP#2 


Figure  5.17  Flight  Point  2  Nonlinear  vs.  Linear  Time  Responses 


Sideslip  Command 

The  aircraft  was  commanded  to  a  5-degree  sideslip  command  for  both  flight  points. 
The  nonlinear  simulation  for  Right  point  1  was  a  little  slow  in  tracking  the  command. 
This  could  have  been  expected  due  to  the  fact  that  the  linear  model  was  also  a  little 
slow  tracking  the  sideslip  command.  The  results  however  were  considered  satisfac¬ 
tory  due  to  the  fact  that  the  simulation  had  no  overshoots  and  the  initial  error  between 
the  command  and  the  aircraft  response  was  not  large.  As  for  the  control  surface  de¬ 
flections,  both  the  linear  and  the  nonlinear  curves  were  smooth  and  not  excessive.  As 
a  matter  of  fact  the  nonlinear  and  the  linear  control  deflections  almost  overlap  com¬ 
pletely. 

The  nonlinear  simulation  for  flight  point  2  was  very  successful  in  tracking  the 
sideslip  command.  It  quickly  converged  to  the  command  without  overshoots.  The 
control  deflection  were  smooth  and  within  tolerances.  The  deflections  were  almost 
on  top  of  one  another  for  the  linear  and  nonlinear  case.  Overall  the  nonlinear  case  was 
successful  for  the  sideslip  command. 

A  graphical  description  of  the  sideslip  command  time  responses  are  found  in  Fig¬ 
ures  5.18-5.21.  In  these  graphs  the  nonlinear  solution  is  always  the  thickest  line.  The 
linear  solution  is  the  thin  line  in  all  cases  except  where  the  input  command  is  also 
plotted.  In  this  case  the  thin  line  is  the  command  input  and  the  dashed  line  is  the  linear 
solution. 


f\  A 


Nonlinear  vs  Linear  Time  Response  for  FP#1 


Betacmd  =  5  degs 
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Figure  5.18  Flight  Point  1  Nonlinear  vs.  Linear  Time  Responses 


Nonlinear  vs  Linear  Time  Response  for  FP#1 


Figure  5.19  Flight  Point  1  Nonlinear  vs.  Linear  Time  Responses 


roll  rate  (deg/sec) 


Nonlinear  vs  Linear  Tinne  Response  for  FP#2  Betacmd  =  5  degs 


time  (sec)  time  (sec) 


Figure  5.20  Right  Point  2  Nonlinear  vs.  Linear  Time  Responses 


Nonlinear  vs  Linear  Time  Response  for  FP#2 


Figure  5.21  Flight  Point  2  Nonlinear  vs.  Linear  Time  Responses 


Chapter  6  Conclusions 

6.1  Summary 

The  objective  of  this  report  was  to  develop  a  total  energy  control  system  (TECS)  for 
the  lateral  dynamics  and  validate  the  design  in  a  nonlinear  simulation  of  the  McDon¬ 
nell  Douglas  F-15.  This  was  accomplished  by  first  analyzing  the  open-loop  charac¬ 
teristics  of  both  the  linear  and  nonlinear  simulations  in  order  to  validate  the  linear 
system.  Once  the  linear  system  was  proven  to  accurately  describe  the  aircraft’s  non¬ 
linear  dynamics,  the  TECS  controller  was  design  for  the  linear  system  through  the 
use  of  the  SANDY  program.  SANDY  is  a  numerical  optimization  program  that  came 
up  with  the  optimum  gains  to  be  used  to  solve  the  control  problem.  These  gains  were 
effective  in  controlling  the  linear  model  both  in  performance  and  robustness.  Next 
these  gains  were  tested  in  the  SIMULE'JK  model  of  the  linear  model  to  make  sure 
there  was  no  extra-problems  shifting  from  the  MATLAB  environment  to  the  SEM- 
ULINK  environment.  This  was  accomplished  locating  several  bugs  in  the  system  be¬ 
fore  the  application  to  the  nonlinear  model  was  accomplished.  The  filial  step  was  to 
test  the  gains  on  the  nonlinear  simulation.  The  gains  were  as  effective  in  tracking  the 
heading  and  sideslip  commands  ^  the  linear  system.  The  end  result  is  an  effective 
Total  Energy  Control  System  for  the  lateral  axis  of  the  McDonnell  Douglas  F-15 
fighter  aircraft. 

6.2  Recommendations  for  Future  Study 

There  are  several  areas  for  future  study.  The  following  is  a  list  where  study  can  be 
furthered. 

1)  Lt.  Dutton  did  the  longitudinal  axis  and  I  did  the  lateral  axis  at  two  specific 
flight  points  [10].  Better  understanding  and  learning  could  be  accomphshed  by  ex¬ 
panding  the  flight  envelope  for  more  than  two  flight  points  and  then  connecting  the 
total  design  with  a  gain  schedule. 

2)  This  thesis  dealt  with  the  McDonnell  Douglas  F-15,  another  future  study 
would  be  to  apply  the  nonlinear  model  framework  to  another  aircraft  using  its  specif¬ 
ic  nonlinear  aerodynamics  and  propulsion  characteristics. 


3)  Control  surface  vibrations  were  discovered  when  applying  the  TECS  control- 
law  to  the  nonlinear  simulation.  Finding  a  way  to  eliminate  these  vibration  is  worth 
investigation. 

These  are  just  a  few  of  the  areas  that  can  be  investigated.  I  hope  this  thesis  pro¬ 
vokes  thought  into  other  aspects  concerning  control  theory.  If  it  does  then  this  paper 
was  a  success  not  only  in  determining  an  effective  control-law  for  the  F-15,  but  also 
it  has  kept  dreams  of  advanced  flight  alive. 
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Appendix  A  F-15  Nonlinear  Simulation  S-Functions 


A.l  S-Funtion  for  Open  Loop  F-15  Model 

The  following  is  a  listing  for  the  MATLAB  code  f25sfn.m  which  is  used  for  the 
open-loop  nonlinear  F-15  simulation. 

%%%  NOmNEAR  F-15  SIMULATION  S-FUNCTION  %%% 
function  [sys,xO]  =  f25sfn(t,x,u,flag) 
global  A  lA 

global  V  alpha  q  theta  p  phi  r  psi  beta  xp  yp  h 

global  DH  DD  DA  DR  PLAPL  PLAPR  PLASYM  alpdot 

global  XX  utrim 

%%%  UPDATE  A,IA-ARRAYS  %%% 

A(829)  =  V;  A(914)  =  alpha;  A(862)  =  q;  A(713)  =  h; 

A(943)  =  theta;  A(861)  =  p;  A(942)  =  phi;  A(715)  =  xp; 

A(863)  =  r;  A(944)  =  psi;  A(915)  =  beta;  A(716)  =  yp; 

A(1402)  =  DH;  A(1403)  =  DD;  A(1401)  =  DA; 

A(1404)  =  DR;  A(1416)  =  PLAPL;  A(1417)  =  PLAPR; 

A(1418)  =  PLASYM; 

%%%  STANDARD  ATMOSPHERE  %%% 

[AMCH,RHO,QBAR,G]  =  atmos(A); 

%%%  A-ARRAY  UPDATE  %%% 

A(825)  =  AMCH;  A(670)  =  RHO;  A(669)  =  QBAR;  A(772)  =  G; 

%%%  STABILITY  AXIS  FORCES  AND  MOMENTS  %%% 
[CLFT,CD,CY,CL,CM,CN,FAX,FAY,FAZ,ALM,AMM,ANM]  =  f25aero(A,- 
lA); 


%%  Transfer  to  the  array  A  (may  not  need  to) 


L  =  QBAR*S*CLFT; 
D  =  QBAR*S*CD; 

Y  =  QBAR*S*CY; 
SumL  =  ALM; 

SumM  =  AMM; 
SumN  =  ANM; 

XT  =  FPX; 

YT  =  FPY; 

ZT  =  FPZ; 


%%%  ANGLE  CALCULATIONS  %%% 

COSTHETA  =  cos(theta); 

SINTHETA  =  sin(theta); 

TANTHETA  =  tan(theta); 

SECTHETA  =  l/cos(theta); 

COSBETA  =  cos(beta); 

SINBETA  =  sin(beta); 

TANBETA  =  tan  (beta); 

COSALPHA  =  cos(alpha); 

SINALPHA  =  sin(alpha); 

COSPHI  =  cos(phi); 

SINPHI  =  sin(phi); 

COSPSI  =  cos(psi); 

SINPSI  =  sin(psi); 

%%%  EQUATIONS  OF  MOTION  %%% 

VI  =  -D*COSBETA+Y*SINBETA+XT*COSALPHA*COSBETA; 
V2  =  YT*SINBETA+ZT*SINALPHA*COSBETA; 

V3  =  -m*g*(COSALPHA*COSBETA*SINTHETA); 

V4  =  -m*g*(-SINBETA*SINPHI*COSTHETA); 
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V5  =  -m*g*(-SINALPHA*COSBETA*COSPHI*COSTHETA); 
all  =  -L+ZT*COSALPHA-XT*SINALPHA; 

al2  =  m*g*(COSALPHA*COSPHI*COSTHETA+SINALPHA*SINTHETA); 
al3  =  q-TANBETA*(p*COSALPHA+r*SINALPHA); 


ql  =  SumL*I2+SumM*I4+SumN*I5-p''2*(Ixz*I4-Ixy*I5); 
q2  =  p*q*(Ixz*I2-Iyz*I4-Dz*I5)-p*r*(Ixy*I2+Dy*I4-Iyz*I5); 
q3  =  q^2*(Iyz*I2-Ixy*I5)-q*r*(Dx*I2-Ixy*I4+Ixz*I5); 
q4  =  -r^2*(Iyz*I2-Ixz*I4); 


pi  =  SumL*Il+SumM*I2+SumN*I3-p''2*(Ixz*I2-Ixy*I3): 
p2  =  p*q*(Ixz*Il-Iyz*I2-Dz*I3)-p*r*(Ixy*Il+Dy*I2-Iyz*I3); 
p3  =  q'^2*(Iyz*Il-Ixy*I3)-q*r*(Dx*Il-Ixy*I2+Ixz*I3); 
p4  =  -r^2*(Iyz*Il-Ixz*I2); 


rl  =  SumL*I3+SumM*I5+SumN*I6-p^2*(Ixz*I5-Ixy*I6); 
r2  =  p*q*(Ixz*I3-Iyz*I5-Dz*I6)-p*r*(Ixy*I3+Dy*I5-Iyz*I6); 
r3  =  q^2*(Iyz*I3-Ixy*I6)-q*r*(Dx*I3-Ixy*I5+Ixz*I6); 
r4  =  -r^2*(Iyz*I3-Ixz*I5); 


bel  =  D=fSINBETA+Y*COSBETA-XT*COSALPHA*SINBETA; 

be2  =  YT*COSBETA-ZT*SINALPHA*SINBETA; 

be3  =  m*g*(COSALPHA*SINBETA*SINTHETA); 

be4  =  m*g*(COSBETA*SINPHI*COSTHETA); 

be5  =  m*g*(-SINALPHA*SINBETA*COSPHPCOSTHETA); 

be6  =  p*SINALPHA-r*COSALPHA; 

xpl  =  COSALPHA*COSBETA*COSTHETA*COSPSI; 

xp2  =  SINBETA*(SINPHI*SINTHETA*COSPSI-COSPHI*SINPSI); 


xp3  =  SINALPHA*COSBETA*(COSPHI*SINTHETA*SINPSI-SINPHI* 
COSPSI); 

ypl  =  COSALPHA*COSBETA*COSTHETA*SINPSI; 
yp2  =  SINBETA*(COSPHI*COSPSI+SINPHI*SINTHETA*SINPSI); 
yp3  =  SINALPHA*COSBETA*(COSPHI*SINTHETA*SINPSI-SINPHI* 
COSPSI); 

hi  =  COSALPHA*COSBETA*SINTHETA; 

h2  =  SINBETA*SINPHI*COSTHETA; 

h3  =  SINALPHA*COSBETA*COSPHI*COSTHETA; 


if  flag  ==  0 


%%%  SYSTEM  CHARACTERISUCS/INITIAL  CONDITIONS  %%% 
sys  =  [12  0  13  8  12  0]; 
xO  =  xx; 

%%%  LONGITUDINAL  STATES  %%% 

V  =  x0(l);  %%%  TOTAL  VEHICLE  VELOCITY  (FT/S)  %%% 
alpha  =  x0(2);  %%%  ANGLE  OF  ATTACK  (RAD)  %%% 
q  =  x0(3);  %%%  PITCH  RATE  (RAD/S)  %%% 
theta  =  x0(4);  %%%  PITCH  ANGLE  (RAD)  %%% 


%%%  LATERAL/DIRECTIONAL  STATES  %%% 
p  =  x0(5);  %%%  ROLL  RATE  (RAD/S)  %%% 
phi  =  x0(6);  %%%  ROLL  ANGLE  (RAD)  %%% 
r  =  x0(7);  %%%  YAW  RATE  (RAD/S)  %%% 
psi  =  x0(8);  %%%  YAW  ANGLE  (RAD)  %%% 
beta  =  x0(9);  %%%  SIDESLIP  ANGLE  (RAD)  %%% 


%%%  EARTH-RELATIVE  POSITION  STATES  %%% 
xp  =  xO(lO);  %%%  X-DIRECnON  POSITION  (FT)  %%% 
yp  =  xO(l  1);  %%%  Y-DIRECTTON  POSITION  (FT)  %%% 
h  =  x0(12);  %%%  ALTITUDE  (FT)  %%% 


elseif  abs(flag)  ==  1 

%%%  LONGITUDINAL  STATES  %%% 

V  =  x(l);  %%%  TOTAL  VEHICLE  VELOCITY  (FT/S)  %%% 
alpha  =  x(2);  %%%  ANGLE  OF  ATTACK  (RAD)  %%% 
q  =  x(3);  %%%  PITCH  RATE  (RAD/S)  %%% 
theta  =  x(4);  %%%  PITCH  ANGLE  (RAD)  %%% 


%%%  LATERAL/DIRECTIONAL  STATES  %%% 
p  =  x(5);  %%%  ROLL  RATE  (RAD/S)  %%% 
phi  =  x(6):  %%%  ROLL  ANGLE  (RAD)  %%% 
r  =  x(7);  %%%  YAW  RATE  (RAD/S)  %%% 
psi  =  x(8);  %%%  YAW  ANGLE  (RAD)  %%% 
beta  =  x(9);  %%%  SIDESLIP  ANGLE  (RAD)  %%% 

%%%  EARTH-RELATIVE  POSITION  STATES  %%% 
xp  =  x(10);  %%%  X-DIRECnON  POSITION  (FT)  %%% 
yp  =  x(ll);  %%%  Y-DIRECnON  POSITION  (FT)  %%% 
h  =  x(12);  %%%  ALTITUDE  (FT)  %%% 


%%%  INPUTS  (U)  %%% 

DH  =  u(l)+utrim(l);  %%%  SYMETRIC  STABILATOR  (DEG)  %%% 
DD  =  u(2);  %%%  DIFFERENTIAL  STABILATOR  (DEG)  %%% 

DA  =  u(3)+utrim(5);  %%%  AILERON  DEFLECTION  (DEG)  %%% 
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DR  =  u(4)+utrim(6);  %%%  RUDDER  DEFLECTION  (DEG)  %%% 
PLAPL  =  u(5)+utrim(2);  %%%  LEFT  PLA  (DEG)  %%% 

PLAPR  =  u(6)+utrim(3);  %%%  RIGHT  PLA  (DEG)  %%% 
PLASYM  =  u(7);  %%%  SYMMETRIC  PLA  (DEG)  %%% 
alpdot  =  u(8);  %%%  AOA  RATE  (RAD/S)  %%% 

%%%  STATE  DERIVATIVES  (dX/dT)  %%% 
sys(l,l)  =  (Vl+V2+V3+V4+V5)/m; 
sys(2,l)  =  al3+(all+al2)/(V*m*COSBETA); 
sys(3,l)  =  (ql+q2+q3+q4)/detr, 
sys(4,l)  =q*COSPHI-r*SINPHI; 
sys(5,l)  =  (pl+p2+p3+p4)/detl; 

sys(6,l)  =  p+q*SINPHPTANTHETA+r*COSPHI*TANTHETA; 
sys(7,l)  =  (rl+r2+r3+r4)/detl; 

sys(8,l)  =  r*COSPHI*SECTHETA+q*SINPHI*SECTHETA; 
sys(9,l)  =  be6+(bel+be2+be3+be4+be5)/(m*V); 
sys(10,l)  =  V*(xpl+xp2+xp3); 
sys(ll,l)  =  V*(ypl+yp2+yp3); 
sys(12,l)  =  V*(hl-h2-h3); 

elseif  flag  ==  3 

%%%  SYSTEM  OUTPUTS  (Y)  %%% 

sys(l,l)  =  x(l); 

sys(2,l)  =  x(2); 

sys(3,l)  =  x(3); 

sys(4,l)  =  x(4); 

sys(5,l)  =  x(5); 

sys(6,l)  =  x(6); 

sys(7,l)  =  x(7); 
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sys(8,l)  =  x(8); 
sys(9,l)  =  x(9); 

sys(10,l)  =  al3+(all+al2)/(V*m*COSBETA); 
sys(ll,l)  =  x(10); 

^(12,l)  =  x(ll); 
sys(13,l)  =  x(12);  ^ - 

else 

%%%  ALL  OTHER  FLAGS  UNDECLARED  %%% 
sys  =  []; 

end 

end 

A.2  S-Fimction  for  Ciosed>Loop  F-15  Model 

The  following  is  a  listing  of  the  MATLAB  s-function  fZSsfncl.m  whicih  is  used 
for  the  closed-loop  nonlinear  F-15  model. 

%%%  NONLINEAR  F-15  SIMULATION  S-FUNCTION  %%% 
function  [sys,xO]  =  f25sfncl(t,x,u,flag) 
global  A  lA 

global  V  alpha  q  theta  p  phi  r  psi  beta  xp  yp  h 

global  DH  PLAPL  PLAPR  FLAPS  DA  DR  PLASYM  DD 

global  XX  utrim 

%%%  UPDATE  A,IA-ARRAYS  %%% 

A(829)  =  V;  A(914)  =  alpha;  A(862)  =  q;  A(713)  =  h; 

A(943)  =  theta;  A(861)  =  p;  A(942)  =  phi;  A(715)  =  xp; 

A(863)  =  r;  A(944)  =  psi;  A(915)  =  beta;  A(716)  =  yp; 
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A(1402)  =  DH;  A(1403)  =  DD;  A(1401)  =  DA; 

A(1404)  =  DR;  A(1416)  =  PLAPL;  A(1417)  =  PLAPR; 
A(1418)  =  PLASYM; 

%%%  STANDARD  ATMOSPHERE  %%% 
[AMCH,RHO,QBAR,G]  =  atmos(A); 


%%%  A-ARRAY  UPDATE  %%% 

A(825)  =  AMCH;  A(670)  =  RHO;  A(669)  =  QBAR;  A(772)  =  G; 

%%%  STABILITY  AXIS  FORCES  AND  MOMENTS  %%% 
[CLFT,CD,CY,CL,CM,CN,FAX,FAY,FAZ,ALM,AMM,ANM]  =  f25aero(A,- 


%%  Transfer  to  the  array  A  (may  not  need  to) 
A(1410)  =  CLFT;  A(14I1)  =  CD;  A(I412)  =  CY; 
A(1413)  =  CL;  A(1414)  =  CM;  A(I415)  =  CN; 
A(748)  =  FAX;  A(749)  =  FAY;  A(750)  =  FAZ; 
A(733)  =  ALM;  A(734)  =  AMM;  A(735)  =  ANM; 


%%%  CALCULATE  PROPULSION  FORCES  %%% 
[FPX,FPY,FPZ,DCL,DCM,DCN,TAUL,TAUR,PLAL,PLAR,COUTlC,COUT 
2C,FIRST]  =  f25eng(A,IA); 


%%  Transfer  to  the  array  A  and  I A 
A(751)  =  FPX;  A(752)  =  FPY;  A(753)  =  FPZ; 
A(736)  =  DCL;  A(737)  =  DCM;  A(738)  =  DCN; 
A(1419)  =  TAUL;  A(1420)  =  TAUR; 

A(1431)  =  PLAL;  A(I432)  =  PLAR; 


A(667)  =  COUTIC;  A(668)  =  COUT2C; 
IA(502)=nRST; 

%%%  A, lA- ARRAY  VAR  NAMES  %%% 

S  =  A(659);  QBAR  =  A(669);  g  =  A(772); 

CLFT  =  A(1410);  CD  =  A(141 1);  CY  =  A(1412); 

ALM  =  A(733);  AMM  =  A(734);  ANM  =  A(735); 

FPX  =  A(751);  FPY  =  A(752);  FPZ  =  A(753); 

Ix  =  A(634);  ly  =  A(635);  Iz  =  A(636); 

Ixz  =  A(637);  Ixy  =  A(638);  lyz  =  A(639); 

II  =  A(632);  12  =  A(631);  13  =  A(630); 

14  =  A(629);  15  =  A(628);  16  =  A(627); 

Dx  =  A(626);  Dy  =  A(625);  Dz  =  A(624); 
detl  =  A(633);  m  =  A(658)/G; 

%%%  CALCULATE  UPDATED  AERO  VALUES  %%% 
L  =  QBAR*S*CLFT; 

D  =  QBAR*S*CD; 

Y  =  QBAR*S*CY; 

SumL  =  ALM; 

SumM  =  AMM; 

SumN  =  ANM; 

XT  =  FPX; 

YT  =  FPY; 

ZT  =  FPZ; 


%%%  ANGLE  CALCULATIONS  %%% 
COSTHETA  =  cos(theta); 

SINTHETA  =  sm(theta); 

TANTHETA  =  tan(theta); 


SECTHETA  =  l/cos(theta); 
COSBETA  =  cos(beta); 
SINBETA  =  sin(beta); 
TANBETA  =  tan(beta); 
COSALPHA  =  cos(alpha); 
SINALPHA  =  sin(alpha); 
COSPHI  =  cos(phi); 
SINPHI  =  sin(phi); 

COSPSI  =  cos(psi); 

SINPSI  =  sin(psi): 


%%%  EQUATIONS  OF  MOTION  %%% 

VI  =  -D*COSBETA+Y*SINBETA+XT*COSALPHA*COSBETA; 

V2  =  YT*SINBETA+ZT*SINALPHA*COSBETA; 

V3  =  -ra*g*(COSALPHA*COSBETA*SINTHETA): 

V4  =  -m*g*(-SINBETA*SINPHI*COSTHETA); 

V5  =  -m*g*(-SINALPHA*COSBETA*COSPHPCOSTHETA); 

all  =  -L+ZT*COSALPHA-XT*SINALPHA; 

al2  =  m*g*(COSALPHA*COSPHI*COSTHETA+SINALPHA*SINTHETA); 
al3  =  q-TANBETA*(p*COSALPHA+r*SINALPHA); 


ql  =  SumL*I2+SumM*I4+SuraN*I5-p^2*(Ixz*I4-Ixy*I5); 
q2  =  p*q*(Ixz*I2-Iyz*I4-Dz*I5)-p*r*(Ixy*I2+Dy*I4-Iyz*I5); 
q3  =  q''2*(Iyz*I2-Ixy*I5)-q*r*(Dx*I2-Ixy*I4+Ixz*I5); 
q4  =  -P2*(Iyz*I2-Ixz*I4); 


pi  =  SumL*Il+SumM*I2+SuraN*I3-p^2*(Ixz*I2-Ixy*I3); 
p2  =  p*q*(Ixz*Il-Iyz*I2-Dz*I3)-p*r*(Ixy*Il+Dy*I2-Iyz*I3); 
p3  =  q^2*(Iyz*Il-Ixy*I3)-q*r*(Dx*Il-Ixy*I2+Ixz*I3); 


p4  =  -r''2*(Iyz*Il-Ixz*I2); 


rl  =  SumL*I3+SuraM*I5+SumN*I6-p'^2*(Ixz*I5-Ixy*I6); 
r2  =  p*q*(Ixz*I3-Iyz*I5-Dz*I6)-p*r*(Ixy*I3+Dy*I5-Iyz*I6); 
r3  =  q''2*(Iyz*I3-Ixy*I6)-q*r*(Dx*I3-Ixy*I5+Ixz*I6); 
r4  =  -r'^2*(Iyz*I3-Ixz*I5); 


bel  =  D*SINBETA+Y*COSBETA-XT*COSALPHA*SINBETA; 

be2  =  YT*COSBETA-ZT*SINALPHA*SINBETA; 

be3  =  m*g*(COSALPHA*SINBETA*SINTHETA); 

be4  =  m*g*(COSBETA*SINPHI*COSTHETA); 

be5  =  m*g*(-SINALPHA*SINBETA*COSPHI*COSTHETA); 

be6  =  p*SINALPHA-r*COSALPHA; 

xpl  =  COSALPHA*COSBETA*COSTHETA*COSPSI; 
xp2  =  SINBETA*(SINPHI*SINTHETA*COSPSI-COSPHI*SINPSI); 
xp3  =  SINALPHA*COSBETA*(COSPHI*SINTHETA*SINPSI-SINPHP 
COSPSI); 

ypl  =  COSALPHA*COSBETA*COSTHETA*SINPSI; 
yp2  =  SINBETA*(COSPHI*COSPSI+SINPHI*SINTHETA*SINPSI); 
yp3  =  SINALPHA*COSBETA*(COSPHI*SINTHETA*SINPSI-SINPHP 
COSPSI); 

hi  =  COSALPHA*COSBETA*SINTHETA; 

h2  =  SINBETA*SINPHPCOSTHETA; 

h3  =  SINALPHA*COSBETA*COSPHPCOSTHETA; 


%disp( ‘before  flag  check’) 
%flag,t,x,u,sys 


%pause 


if  flag  ==  0 


%%%  SYSTEM  CHARACTERISTICS/INITIAL  CONDITIONS  %%% 
sys  =  [12  0  14  6  12  0]; 
xO  =  xx; 

%%%  LONGITUDINAL  STATES  %%% 

V  =  x0(l);  %%%  TOTAL  VEHICLE  VELOCITY  (FT/S)  %%% 
alpha  =  x0(2);  %%%  ANGLE  OF  ATTACK  (RAD)  %%% 
q  =  x0(3);  %%%  PITCH  RATE  (RAD/S)  %%% 
theta  =  x0(4);  %%%  PITCH  ANGLE  (RAD)  %%% 

%%%  LATERAL/DIRECTIONAL  STATES  %%% 
p  =  x0(5);  %%%  ROLL  RATE  (RAD/S)  %%% 
phi  =  x0(6);  %%%  ROLL  ANGLE  (RAD)  %%% 
r  =  x0(7);  %%%  YAW  RATE  (RAD/S)  %%% 
psi  =  x0(8);  %%%  YAW  ANGLE  (RAD)  %%% 
beta  =  x0(9);  %%%  SIDESLIP  ANGLE  (RAD)  %%% 

%%%  EARTH-RELATIVE  POSITION  STATES  %%% 
xp  =  xO(lO);  %%%  X-DIRECnON  POSITION  (FT)  %%% 
yp  =  x0(l  1);  %%%  Y-DIRECTTON  POSITION  (FT)  %%% 
h  =  x0(12);  %%%  ALTITUDE  (FT)  %%% 

elseif  abs(flag)  ==  1 


%%%  SYSTEM  STATES  (X)  %%% 
%%%  LONGITUDINAL  STATES  %%% 


V  =  x(l);  %%%  TOTAL  VEHICLE  VELOCITY  (FT/S)  %%% 
alpha  =  x(2);  %%%  ANGLE  OF  ATTACK  (RAD)  %%% 
q  =  x(3);  %%%  PITCH  RATE  (RAD/S)  %%% 
theta  =  x(4);  %%%  PITCH  ANGLE  (RAD)  %%% 

%%%  LATERAL/DIRECTIONAL  STATES  %%% 
p  =  x(5);  %%%  ROLL  RATE  (RAD/S)  %%% 
phi  =  x(6);  %%%  ROLL  ANGLE  (RAD)  %%% 
r  =  x(7):  %%%  YAW  RATE  (RAD/S)  %%% 
psi  =  x(8);  %%%  YAW  ANGLE  (RAD)  %%% 
beta  =  x(9);  %%%  SIDESLIP  ANGLE  (RAD)  %%% 

%%%  EARTH-RELATIVE  POSITION  STATES  %%% 
xp  =  x(10);  %%%  X-DIRECTION  POSITION  (FT)  %%% 
yp  =  x(l  1);  %%%  Y-DIRECnON  POSITION  (FT)  %%% 
h  =  x(12);  %%%  ALTITUDE  (FT)  %%% 


%%%  INPUTS  (U)  %%% 

DH  =  u(l)+utrim(l);  %%%  SYMETRIC  STABILATOR  DEFLECTION  (DEG) 

%%% 

PLAPL  =  u(2)+utrim(2);  %%%  LEFT  ENGINE  POWER  LEVEL  ANGLE 
(DEG)  %%% 

PLAPR  =  u(3)+utrim(3);  %%%  RIGHT  ENGINE  POWER  LEVEL  ANGLE 
(DEG)  %%% 

FLAPS  =  u(4)-mtrim(4);  %%%  Flaps  (DEG)  %%% 

DA  =  u(5)+utrim(5);  %%%  DIFFERENTIAL  AILERON  DEFLECTION 
(DEG)  %%% 

DR  =  u(6)+utrim(6);  %%%  RUDDER  DEFLECTION  (DEG)  %%% 

PLASYM  =  0; 


%%%  STATE  DERIVATIVES  (dX/dT)  %%% 
sys(l,l)  =  (Vl+V2+V3+V4+V5)/m; 
sys(2il^^  al3+(all+al2)/(V*m*COSBETA); 
sys(3,l)  =  (ql+q2+q3+q4)/detl; 
sys(4,l)  =  q*COSPHI-r*SINPHI; 
sys(5,l)  =  (pl+p2+p3+p4)/detl; 

sys(6,l)  =  p+q*SINPHI*TANTHETA+r*COSPHPTANTHETA; 
sys(7,l)  =  (rl+r2+r3+r4)/detl; 

sys(8,l)  =  r*COSPHPSECTHETA+q*SINPHI*SECTHETA; 
sys(9,l)  =  be6+(bel+be2+be3+be4+be5)/(m*V); 
sys(10,l)  =  V*(xpl+xp2+xp3); 
sys(ll,l)  =  V*(ypl+yp2+yp3); 
sys(12,l)  =  V*(hl-h2-h3); 

elseif  flag  ==  3 

%%%  SYSTEM  OUTPUTS  (Y)  %%% 

sys(l,l)  =  (Vl+V2+V3+V4+V5)/m/g;  %%%  V_DOT/G  %%% 

sys(2,l)  =  x(l)-xx(l);  %%%  V  %%% 

sys(3,l)  =  x(12)-xx(12);  %%%  H  %%% 

sys(4,l)  =  x(4)-xx(4)-x(2)+xx(2);  %%%  GAMMA  %%% 

sys(5,l)  =  x(3)-xx(3);  %%%  Q  %%% 

sys(6,l)  =  x(4)-xx(4);  %%%  THETA  %%% 

sys(7,l)  =  x(2)-xx(2);  %%%  ALPHA  %%% 

sys(8,l)  =  x(5)-xx(5);  %%%  P  %%% 

sys(9,l)  =  x(6)-xx(6);  %%%  PHI  %%% 

ifV<  le-3 


79 


sys(10,l)  =  0; 
else 

sys(10,l)  =  be6+(bel+be2+be3+be4+be5)/(m*V);  %%%  BETA_DOT  %%% 
end 

sys(ll,l)  =  x(9)-xx(9);  %%%  BETA  %%% 
sys(12,l)  =  x(8)-xx(8);  %%%  PSI  %%% 

sys(13,l)  =  r*COSPHI*SECTHETA+q*SINPHPSECTHETA;  %%%  PSI_- 
DOT  %%% 


sys(14,l)  =  x(7)-xx(7);  %%%  R  %%% 


else 


%%%  ALL  OTHER  FLAGS  UNDECLARED  %%% 
sys  =  []; 


end 

end 


# 


Appendix  B  Linearized  State-Space  Models 

The  following  are  state-space  models  of  the  linearized  F-15  model  for  the  two  flight 
points. 

States  =  [AV(fl/sec),  Aa(rad)  ,Aq(rad/sec),  A0(rad),  Ap(rad/sec),  A(p(rad),  Ar(rad/ 
sec),  A\|r(rad),  AP(rad),  Ah(ft)]^ 

Inputs  =  [ASnCdeg),  A5pLAL(deg),  A5pLAR(deg),  ASpiapCdeg),  A5a(deg),  A5i.(deg)]'^ 


Outputs  =  [AV/g,  AV(ft/sec),  Dh(ft),  AY(rad),  Aq(rad/sec),  A6(rad),  Aa(rad), 
Ap(rad/sec),  A(p(rad),  AP (rad/sec),  AP(rad),  A\j/(rad),  Aijt  (rad/sec),  Ar(rad/sec)]^ 

%7c%  F-15  State  Space  Model  %%% 

%%%  FLIGHT  POINT  #2  %%% 

%%%  9,800  ft,  0.5  Mach  %%% 


al  =  [-1.3658e-02  -3.4605e+00  0  -3.2144e+01  0  -1.2813e-05 

-2.1654e-04  -7.7897e-01  l.OOOOe+00  -2.9875e-05  0  -2.9715e-07 

4.0744e-04  -5.8027e+00  -2.50I3e+00  -6.5003e-04  3.1496e-08  4.0402e-08 

0  0  l.OOOOe+00  0  0  3.2285e-20 

7.4534e-22  2.4167e-19  -1.5331e-19  0  -2.2133e+00  0 

0  0  -8.2426e-20  -3.2493e-20  l.OOOOe+00  -2.2495e-20 

-9.4443e-23  9.4468e-20  1.8812e-20  0  -7.4308e-02  0 

0  0  -1.0274e-18  -2.6070e-21  0  -2.8039e-19 

1.1581e-22  -2.6938e-20  0  4.3973e-2I  7.9981e-02  5.9435e-02 

2.4600e-04  -5.3908e+02  0  5.3908e+02  0  2.1489e-04]; 


a2  =  [0  0  -1.0357e-03  -8.7382e-05 


0  0  1.6619e-12  1.3835e-08 
-3.1496e-08  0  -4.5464e-13  -3.7846e-09 
1.0241e-18  0  00 
1.3971e+00  0-2.7048e+01  0 
8.0486e-02  0  0  0 
-5.7386e-01  0  4.6745e+00  0 
1.0032e+00  0  0  0 

-9.9680e-01  0  -1.9213e-01  -1.3645e-27 
0  0  -6.6284e-07  0]; 


a  =  [ala2]; 

b  =  [-7.1779e-02  1.0670e-01  1.0147e-01  0  0  0 
-L5698e-03  -1.5881e-05  -1.5103e-05  0  0  0 
-1.5079e-01  4.3444e-06  4.1314e-06  0  0  0 
000000 

0  0  0  0  1.6930e-01  2.6623e-02 
000000 

0  0  0  0  2.2413e-03  -4.8714e-02 
000000 

-1.1210e-24  1.6663e-24  L5846e-24  0 -3.8274e-05  6.3187e-04 
0  0  0  0  0  0]; 


cl  =  [-4.2492e-04  -1.0766e-01  0  -l.OOOOe+00  0  -3.9862e-07 
l.OOOOe+00  0  0  0  0  0 
000000 

0  -l.OOOOe+00  0  l.OOOOe+00  0  0 
0  0  l.OOOOe+OO  0  0  0 
0  0  0  l.OOOOe+00  0  0 
0  l.OOOOe+00  0  0  0]; 
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03  3.15660-03  0  00 
000000 
000000 
000000 
000000 
000000 
000000 
000000 

-L1210e-24  1.6663e-24  1.5846e-24  0 -3.8274e-05  6.3187e-04 
000000 
000000 
000000 
0  0  0  0  0  0]; 


%%%  F-15  State  Space  Model  %%% 

%%%  FLIGHT  POINT  #2  %%% 
%%%  30,000  ft,  0.5  Mach  %%% 


al  =  [-2.0004e-02  -3.8596e+01  0  -3.2083e+01  0  -2.9032e-05 
-2.6078e-04  -3.5479e-01  l.OOOOe+00  -2.8115e-04  0  -3.1179e-07 
2.6521e-04  -2.3011e+00  -1.1476e+00  -6.4020e-04  3.1496e-08  1.1296e-08 
0  0  l.OOOOe+00  0  0  1.9763e-20 
1.3725e-21  -3.2396e-19  -9.3566e-20  0  -1.0499e+00  0 
0  0  -1.0213e-20  -2.0490e-20  l.OOOOe+00  2.6320e-23 
-1.4110e-22  1.7417e-19  5.3361e-20  0  -7.5227e-03  0 
0  0  -5.4237e-20  -3.8585e-21  0  1.3977e-22 
1.1038e-23-7.7720e-20  0-1.1030e-21  1.8428e-01  6.3359e-02 
4.1000e-03  -4.9731e+02  0  4.9731e+02  0  4.5002e-04]; 


a2  =  [0  0  -4.6970e-04  -2.1911e-04 
0  0  5.3732e-12  8.3523e-08 
-3.1496e-08  0  -4.9397e-13  -7.6785e-09 
5.3267e-20  0  0  0 
1.2101e+00  0-1.3163e+010 
L9174e-01  0  0  0 
-3.0149e-01  0  1.4547e+00  0 
1.0182e+00  0  0  0 

-9.8287e-01  0  -9.4448e-02  -1.1953e-26 
0  0 -1.0194e-05  0]; 


a  =  [al  a2]; 


b  =  [-1.4240e-01  3.8878e-02  3.6972e-02  0  0  0 
-7.0387e-04  -1.4657e-05  -1.3939e-05  0  0  0 
-5.5398e-02  1.3475e-06  1.2814e-06  0  0  0 
000000 

0  0  0  0  5.5926e-02  5.4716e-03 
000000 

0  0  0  0  1.4022e-04  -2.1371e-02 
000000 

-7.7682e-24  2.1208e-24  2.0169e-24  0  2.9883e-06  3.1658e-04 
0  0  00  00]; 


cl  =  [-6.2352e-04  -1.2030e+00  0  -9.9999e-01  0  -9.0492e-07 
l.OOOOe+00  00  00  0 
000000 

0  -l.OOOOe+00  0  l.OOOOe+00  0  0 
0  0  l.OOOOe+00  0  0  0 
0  0  0  l.OOOOe+00  0  0 


0  l.OOOOe+00  0  0  0  0 
0  0  0  0  l.OOOOe+00  0 
0  0  0  0  0  l.OOOOe+00 

1.1038e-23  -7.7720e-20  0  -1.1030e-21  1.8428e-01  6.3359e-02 

booooo 

000000 

0  0  -5.4237e-20  -3.8585e-21  0  1.3977e-22 
00  00  00]; 


c2  =  [0  0  -1.4640e-05  -6.8295e-06 
0000 

0  0  0  l.OOOOe+00 
0000 
0000 
0000 
0000 
0000 
0000 


-9.8287e-01  0  -9.4448e-02  -1.1953e-26 
0  0  l.OOOOe+00  0 
0  l.OOOOe+00  0  0 
1.01826+00  00  0 
l.OOOOe+00  0  0  0]; 


c  =  [cl  c2]; 


d  =  [-4.43856-03  1.21186-03  1.15246-03  0  0  0 
000000 
000000 
000000 


000000 
000000 
000000 
00  00  0  0 
00t)000 

-7.7682e-24  2.1208e-24  2.0169e-24  0  2.9883e-06  3.1658e-04 
000000 
000000 
000000 
0  0  0  0  0  0]; 


Appendix  C  Files  used  by  SANDY  to  Acquire  the  Optimum  Gains 
C.1  Flight  Point  1 

This  file  is  zzlatTECSl.dat  which  contains  the  data  used  by  the  SANDY  program  to 
come  up  with  the  optimum  gains  for  the  first  flight  point. 

‘Nfcmax’  11 
2000 

‘Npm’  1  1 
1 

‘Tf  1  1 
100 

‘Tfctor’  1  1 
2 

‘Wp’  1  1 
1 


‘F  9  9 

-0.094448  0.18428  0.063359 
-0.98287  0.00017122959  0.018140034 
0  0.00164792149893461  0.010050294940281 
-13.163-1.0499  0 
1.2101  3.2045598  0.31352268 
0  0.229667019846649  1.4006864337934 
0  1  2.632e-23 
0.19174  0  0 
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1.4547  -0.0075227  0 
-0.30149  0.008034606  -1.2245583 
0  -0.0253814946266748  -0.154795909385342 
000 
0-20  0 
0  -0  -0 
000 
0  0-20 
0-0  -0 
000 
100 
000 
000 
000 
00  1 
000 
000 

0  -0.080656  -0.5396 


‘G’  9  2 
00 
00 
00 
00 
20  0 
0  20 
00 
00 
00 
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‘Gam’  9  3 
000 
000 
000 
000 
000 
000 
000 
000 
100 


‘Hs’  6  9 
000 
100 
000 

-0.094448  0.18428  0.063359 
-0.98287  0.00017122959  0.018140034 
0  0.00164792149893461  0.010050294940281 
0  10 
000 
000 
00  1 
000 
000 
100 
000 

0  -0.0174479237139443  -0.106410881546258 
000 
000 
100 
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‘Dsu’  6  2 
00 
00 
00 
00 
00 
00 

‘Dsw’  6  3 
000 
000 
000 
000 
0-1  0 
00-1 


‘He’  3  9 

-0.00538843614841794  -10.4180780028834  -0.0171376540600306 
-8.37362741591059  4.61237608695652e-05  0.00488634347826087 
0  9.4017022855056e-05  0.000573388240708870 
100 
000 

0  -0.0174479237139443  -0.106410881546258 
000 
000 
100 


‘Dcu’  3  2 


‘Wspec’  3  6 
000 
00  1 
0  0  0.1 
10.10 
0  0  0.1 
10.10 

‘Alpha’  1  1 
0 

‘Sigmamax’  1  1 
-0.15 


‘Zetamin’  1  1 
0.7 


‘Q’  3  3 
00  0 
0  20  0 


‘nDircu’  1  1 
0 


‘Dircu’  0  0 


‘nDircy’  1  1 
0 


‘Dircy’  0  0 


‘A’ 2  2 
00 
00 

‘B’2  6 
-1  -1  0 

0  -3.65667335253509  -1.01743461636475 
1-10 

0  -3.65667335253509  1.01743461636475 


‘C2  2 

-0.253683784709487  0 
0  0.10496436553704 

‘D’  2  6 

1.3514438906174  1.3514438906174  -0.560850786142026 
0.164780307434178  0  0 
-0.371252547475194  -1.82026154119372  0 
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m 


‘nAid’  1  1 
0 

‘nBid’  1  1  ‘‘ 
4 

‘Bid’  4  5 

15  1 

-10000  10000 

16  1 

-10000  10000 

25  1 

-10000  10000 

26  1 

-10000  10000 
‘nCid’  1  1 
2 

‘Cid’  2  5 
1  1  1 

-10000  10000 
22  1 

-10000  10000 
‘nDid’  1  1 
6 

‘Did’  6  5 
1  1  1 

-10000  10000 
12  1 

-10000  10000 
13  1 
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-10000  10000 
14  1 

-10000  10000 
2  1  1 

-10000  10000 
22  1 

-10000  10000 
‘Nclin’  1  1 
3 

‘nLincoef  3  1 
2 
2 
2 

‘Lincoef  6  4 
198  1 
5 

-1  98  2 

5 

198  1 

6 

198  2 
6 

1  100  1 
1 

-1  100  1 
2 

‘Linbnds’  3  2 
00 
00 
00 


‘Nudnlc’  1  1 
0 


C.2  Flight  Point  2 

Here  is  the  MATLAB  file  zzlatTECS2.dat  which  the  SANDY  program  used 
calculate  the  optimum  gains  for  the  second  flight  point. 

‘Nfcmax’  1  1 
2000 

‘Npm’  1  1 
1 

Tf  1  1 
100 

‘Tfctor’  11 
2 

‘Wp’  1  1 
1 

‘F’  9  9 

-0.094448  0.18428  0.063359 
-0.98287  0.00017122959  0.018140034 
0  0.00164792149893461  0.010050294940281 
-13.163-1.0499  0 
1.2101  3.2045598  0.31352268 
0  0.229667019846649  1.4006864337934 
0  1  2.632e-23 
0.19174  0  0 
0-0  -0 

1.4547  -0.0075227  0 

-0.30149  0.008034606  -1.2245583 

0  -0.0253814946266748  -0.154795909385342 
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# 


0  -20  0 
0-0  -0 
000 

0  0  -20  0  -0  -0 
000 
100 
000 
000 
000 
00  1 
000 
000 

0  -0.080656  -0.5396 


‘G’  9  2 
00 
00 
00 
00 
20  0 
0  20 
00 
00 
00 

‘Gam’  9  3 
000 

0  00  000 
000 
00  0 
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000 

000 

000 

100 
‘Hs’  6  9 
000 
100 
00  0 

-0.094448  0.18428  0.063359 
-0.98287  0.00017122959  0.018140034 
0  0.00164792149893461  0.010050294940281 
0  10 
000 
000 
00  1 
000 
000 
100 
000 

0  -0.0174479237139443  -0.106410881546258 
000 
000 
100 

‘DOO 

00 

00 

00 

00 
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‘Dsw’  6  3 
000 
000 

000 

000 

0-1  0 

0  0-1 
‘Hg’  3  9 

-0.00538843614841794  -10.4180780028834  -0.0171376540600306 
-8.37362741591059  4.61237608695652e-05  0.00488634347826087 
0  9.40170228550566-05  0.000573388240708870 

100 

000 

0  -0.0174479237139443  -0.106410881546258 

000 
000 
100 
‘Dcu’  3  2 
00 
00 
00 

‘Dew’  3  3 
000 


‘Wspec’  3  6 
000 
00  1 
0  0  0.1 
10.10 
0  0  0.1 
1  0.1  0 

‘Alpha’  1  1 
0 

‘Sigmamax’  1  1 
-0.15 

‘Zetamin’  1  1 
0.7 

‘Q’  3  3 
000 
0  20  0 
00  1 

‘R’  2  2 
00 
00 

‘nDircu’  1  1 
0 

‘Dircu’  0  0 


‘nDircy’  1  1 
0 

‘Dircy’  0  0 
‘A’  2  2 
00 
00 

‘B’2  6 
-1  -1  0 

0  -3.65667335253509  -1.01743461636475 
1-10 

0  -3.65667335253509  1.01743461636475 


‘C’  2  2 

-0.253683784709487  0 
0  0.10496436553704 


‘D’  2  6 

1.3514438906174  1.3514438906174  -0.560850786142026 
0.164780307434178  0  0 
-0.371252547475194  -1.82026154119372  0 
000 


‘nAid’  1  1 
0 

‘Aid’  0  0 
‘nBid’  1  1 
4 

‘Bid’  4  5 
15  1 

-10000  10000 


-10000  10000 

25  1 

-10000  10000 

26  1 

-10000  10000 
‘nCid’  1  1 
2 

‘Cid’  2  5 
1  1  1 

-10000  10000 
22  1 

-10000  10000 
‘nDid’  1  1 
6 

‘Did’  6  5 
1  1  1 

-10000  10000 
12  1 

-10000  10000 

1  3  1 

-10000  10000 
14  1 

-10000  10000 

2  1  1 

-10000  10000 
22  1 

-10000  10000 
‘Nclin’  1  1 


‘nLincoef  3  1 
2 
2 
2 

‘Lincoef  6  4 

1  98  1 
5 

-1  98  2 

5 

198  1 

6 

198  2 
6 

1  100  1 
1 

-1  100  1 

2 

‘Linbnds’  3  2 
00 
00 
00 

‘Nudnlc’  1  1 


Appendix  D  SANDY  Gains 


D.l  SANDY  Gains  for  Flight  Point  1 

Kp  =  -7.8735 

K(p  = -5.5103 
Kp  = -0.1563 
Ky  =  -0.4172 
K,=  l 
Ka=l 

Kri  =  41.6557 
Krp  =  -95.7255 
Kyp=  1.0616 

Kyi  =  0.6251 
-0.3364 


D.2  SANDY  Gains  for  Flight  Point  2 
Kp  =  1.7888 

K(p  =  0.5208 
Kp  = -2.2818 
Ky  = -0.9177 
Kr=l 
Ka=l 

Kri  =  0.2142 
Krp  =  11.7688 
Kyp=  1.0322 
Kyi  =  0-4382 
-0.5067 


Appendix  E  Operating  Instructions 
E.l  Nonlinear  Open-Loop  Simulation 

The  following  is  a  step-by-step  instruction  needed  to  execute  the  nonlinear  open-loop 
F-15  simulation; 

1)  Start  MATLAB 

2)  Using  an  available  editor,  modify  trimmod.m  in  order  to  chose  the  desired 
flight  point. 

3)  On  the  MATLAB  window  type  fZSload  in  order  to  activate  f251oad.m. 

4)  Type  f25sim  in  MATLAB  in  order  to  bring  up  the  SIMULINK  environment. 

5)  Select  Simulation  menu,  select  Parameters  and  set  the  desired  time  span  as 
well  as  the  method  of  calculation.  This  report  used  the  EULER  method  with  a  mini¬ 
mum  time  step  of  0.001  seconds  in  order  to  avoid  numerical  inconsistencies. 

6)  Chose  the  desired  input  by  clicking  on  the  respective  input  and  typing  in  the 
magnitude  of  the  input. 

7)  Select  the  Simulation  menu  and  choose  Start  to  start  the  simulation. 

8)  Once  the  simulation  is  complete,  the  selected  parameters  of  the  simulation  will 
be  stored  in  the  MATLAB  workspace.  Saving  the  MATLAB  variables  will  make 
them  available  for  future  usage. 

E.2  Nonlinear  Closed-Loop  Simulation 

The  following  is  a  summation  of  the  commands  used  to  run  the  nonlinear  closed-loop 
simulation. 

1)  Start  MATLAB 

2)  Using  an  available  editor,  modify  trimmod.m  in  order  to  choose  the  desired 
flight  point. 

3)  On  the  MATLAB  window  type  f251oad  in  order  to  activate  f251oad.m. 

4)  The  simulations  for  each  flight  point  have  been  saved  under  different  file 
names  for  convenience.  In  order  to  run  Flight  Point  1  type  simu539n  and  for  Flight 
Point  2  type.  simu497n.  These  commands  will  set  up  the  SIMULINK  environment. 

5)  Select  Simulation  menu,  select  Parameters  and  set  the  desired  time  span  as 


well  at  the  method  of  calculation.  This  report  used  the  EULER  method  with  a  mini¬ 
mum  time  step  of  0.001  seconds  in  order  to  avoid  numerical  inconsistencies. 

6)  Choose  the  desired  input  by  clicking  on  the  respective  input  and  typing  in  the 
magnitude  of  the  input. 

7)  Select  the  Simulation  menu  and  choose  Start  to  start  the  simulation. 

8)  Once  the  simulation  is  complete,  the  selected  parameters  of  the  simulation  will 
be  stored  in  the  MATLAB  workspace.  Saving  the  MATLAB  variables  will  make 
them  available  for  future  usage 

E.3  Linear  Closed-Loop  Simulation 

Below  is  a  summary  of  the  necessary  commands  to  run  the  linearized  closed-loop 
simulation. 

1)  Start  MATLAB 

2)  The  simulations  for  each  flight  point  have  been  saved  under  different  file 
names  for  convenience.  In  order  to  run  Flight  Point  1  type  simu539  and  for  Flight 
Point  2  type  simu497.  These  commands  will  set  up  the  SIMULINK  environment. 

3)  Select  Simulation  menu,  select  Parameters  and  set  the  desired  time  span  as 
well  at  the  method  of  calculation.  This  report  used  the  LINSIM  method  with  a  mini¬ 
mum  time  step  of  0.01  seconds  in  order  to  avoid  numerical  inconsistencies. 

4)  Choose  the  desired  input  by  clicking  on  the  respective  input  and  typing  in  the 
magnitude  of  the  input. 

5)  Select  the  Simulation  menu  and  choose  Start  to  start  the  simulation. 

6)  Once  the  simulation  is  complete,  the  selected  parameters  of  the  simulation  will 
be  stored  in  the  MATLAB  workspace.  Saving  the  MATLAB  variables  will  make 
them  available  for  future  usage. 


